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ABSTRACT 


Small  symmetric  dual  reflector  antennas  suffer  from  low 
efficiency  due  to  subreflector  blockage  of  the  main  reflector 
and  subreflector  scattering.  These  can  be  reduced  by  slicing 
the  main  dish  and  translating  its  rotational  axis,  along  with 
modifying  the  subreflector  geometry. 

This  type  of  design  is  usually  applied  to  low-frequency 
reflectors,  but  high-frequency  analysis  techniques  are  used. 
Consequently  the  agreement  between  measured  and  computed  data 
is  not  as  good  as  it  would  be  for  a  rigorous  solution  such  as 
the  method  of  moments. 

This  thesis  modifies  an  existing  method  of  moments  computer 
code  to  handle  the  displaced  axis  geometry,  and  computes  the 
radiation  pattern  and  the  efficiency  of  this  antenna  as  a 
function  of  geometrical  and  electrical  design  parameters. 
Optimum  configurations  are  identified  for  several  feed  types. 
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I .  INTRODUCTION 

A.  DUAL-REFLECTOR  ANTENNA  SYSTEMS 

The  paraboloidal  antenna  with  a  feed  at  the  focus  does  not 
allow  much  control  of  the  power  distribution  over  the  aperture 
surface,  except  for  what  can  be  accomplished  by  changing  the 
focal  length  and  feed  pattern.  The  introduction  of  a  second 
reflecting  surface,  such  as  in  the  Cassegrain  or  Gregorian 
systems,  allows  more  control  of  the  aperture  distribution 
because  of  the  extra  degree  of  freedom  that  a  second 
reflecting  surface  gives.  This  additional  control  over  the 
aperture  field  distribution  is  obtained  by  shaping  both  the 
subreflector  and  the  main  reflector  so  as  to  change  the  power 
distribution  on  the  main  reflector  aperture  but  yet  maintain 
the  required  phase  distribution.  A  single  reflector  does  not 
allow  both  the  power  and  phase  distribution  to  be  varied 
independently  [Ref.  1]. 

In  general  these  types  of  antennas  provide  a  variety  of 
benefits,  such  as  the 

•  ability  to  place  the  feed  in  a  convenient  location 

•  reduction  of  spillover  and  minor  lobe  radiation 

•  ability  to  obtain  an  equivalent  focal  length  much 
greater  than  the  physical  length 

•  capability  for  scanning  and/or  broadening  of  the  beam 
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by  moving  one  of  the  reflecting 
[Ref.  2]. 


surfaces 


The  Cassegrain  antenna,  which  is  shown  in  Figure  1  and 
derived  from  telescope  design  [Ref.  3],  is  the  most 
common  antenna  using  multiple  reflectors.  The  feed  illuminates 
the  hyperboloidal  subreflector,  which  in  turn  illuminates  the 
paraboloidal  main  reflector.  The  feed  is  placed  at  one  focus 
of  the  hyperboloid  and  the  paraboloid  focus  at  the  other.  A 
similar  system  is  the  Gregorian,  which  uses  an  ellipsoidal 
subreflector  in  place  of  the  hyperboloid  as  shown  in  Figure  2. 


Figure  1:  CLASSICAL  CASSEGRAIN 
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Figure  2:  CLASSICAL  GREGORIAN 

Aperture  blocking  can  be  large  for  these  types  of  antennas.  It 
may  be  minimized  by  choosing  the  diameter  of  the  subreflector 
equal  to  that  of  the  feed.  This  "minimum  blocking  condition" 
is  derived  using  geometrical  optics  and  does  not  always  give 
good  results  for  electrically  small  antennas. 

In  the  general  dual-reflector  case,  blockage  can  be 
eliminated  by  either  of  two  ways: 

(1)  offsetting  both  the  feed  and  the  subreflector 
[Ref.  4]  and 

(2)  displacing  the  main  reflector  and  compensating  for 
the  displacement  by  reorienting  the  subreflector. 

In  the  latter  case  rotational  symmetry  is  maintained. 


B.  DISPLACED  AXIS  DUAL-REFLECTOR  ANTENNA  SYSTEMS 

Mobile  satellite  communication  terminals  require  small 
antennas  with  high  aperture  efficiencies.  The  Cassegrain  and 
Gregorian  antenna  types,  which  are  simple  and  therefore 
inexpensive,  are  generally  not  suitable  for  this  purpose  since 
their  subreflectors  have  a  relatively  large  blocking  area, 
causing  a  reduction  in  gain.  In  other  words,  they  have  low 
efficiency  for  the  sizes  of  interest  (  20  to  30  wavelengths) . 

One  method  of  reducing  blockage,  yet  maintaining 
rotational  symmetry,  is  to  displace  the  main  reflector  away 
from  the  axis  of  symmetry  and  then  adjust  the  subreflector  so 
that  a  plane  wave  front  is  achieved  in  the  aperture  plane.  A 
displaced-axis  antenna  with  a  subreflector  that  is  a  portion 
of  a  hyperboloid  will  be  called  an  ATH  for  short,  and  one  with 
elliptical  based  subreflector  an  ATE  (or  ADE) .  Geometrical 
optics  analysis  of  the  ATE  and  ATH  can  be  found  in  [Ref.  5]. 

As  mentioned  earlier,  the  ADE  antenna  may  be  considered  as 
a  special  case  of  the  generalized  Gregorian  system.  However  in 
this  case  the  focal  axis  of  the  main  parabolic  reflector  is 
displaced  from  the  axis  of  symmetry,  which  contains  the  prime 
focus  of  the  elliptical  subreflector  as  shown  in  Figure  3. 
The  locus  of  the  secondary  foci  of  the  subreflector  forms  a 
ring  through  which  all  the  rays  pass,  which  coincides  with  the 
ring  focus  of  the  main  reflector.  The  parameters  of  the 
ellipse  and  the  parabola  are  chosen  so  that  the  central  ray 
from  the  prime  focus  reaches  the  outer  edge  of  the  main 
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reflector,  while  the  inner  ray  clears  the  outer  edge  of  the 
subreflector  after  reflection  from  the  main  reflector. 


The  principal  merits  of  this  geometry  are  that  the  rays  are 
not  reflected  into  the  prime  feed  from  the  subreflector,  nor 
back  into  the  subreflector  from  the  main  reflector.  Also,  the 
aperture  illumination  for  the  radiated  wave  is  more  uniform 
than  that  in  the  standard  Cassegrain  or  Gregorian 
configuration.  All  rays  reflected  from  the  paraboloid  miss  the 
subreflector,  leading  to  a  higher  aperture  efficiency.  Another 
advantage  of  the  ADE  antenna  is  that  the  F/D  ratio  of  the  main 
reflector  can  be  made  smaller  than  that  for  a  regular  design, 
leading  to  a  compact  antenna  design  with  a  reduction  in  far- 
out  sidelobes  [Ref.  6]. 

To  summarize,  the  advantages  of  this  type  of  design  are 
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•  the  main  reflector  may  be  made  as  small  as  ten 

wavelengths  while  retaining  a  relatively  high 

efficiency 

•  a  low  feed  input  reflection  coefficient  can  be 

achieved 

•  small  sidelobe  and  gain  degradation  due  to  aperture 
blockage 

•  elimination  of  the  spares  which  normally  hold  the 

subreflector  (the  subreflector  can  usually  be 

integrated  into  the  feed  design  for  small  antennas) . 

•  subreflectors  which  can  be  several  times  smaller  than 
those  in  Cassegrain  antennas. 

C.  SCOPE  OF  THE  STUDY 

Chapter  II  discusses  the  derivation  and  the  MM  solution  of 
the  electric  field  integral  equation  for  a  general  body  of 
revolution  (BOR)  and  measurements  matrices  are  developed. 
Chapter  III  discusses  the  MM  solution  for  the  ADE  antenna 
geometry  and  evaluates  the  excitation  vector  for  an  ideal 
feed.  Chapter  IV  explains  the  computer  program  that  has  been 
developed  based  on  an  existing  code  written  by  Mautz  and 
Harrington.  Chapter  V  discusses  the  results  obtained  by  the 
program  and  presents  some  guidelines  for  optimizing  the 
geometry.  Chapter  VI  discusses  conclusions  based  on  the 
results  and  presents  some  recommendations  for  further 
research . 
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II. 


INTEGRAL  EQUATIONS  AND  THE  MOMENT  METHOD 


A.  ELECTRIC  FIELD  INTEGRAL  EQUATION 

In  the  following  analysis  time-harmonic  field  quantities 
are  assumed.  Phasor  representations  of  the  fields  are  used 
with  the  e-**1  time  dependence  suppressed. 

The  electric  field  integral  equation  (EFIE)  is  based  on 
the  boundary  condition  that  the  total  tangential  electric 
field  on  a  perfectly  electric  conducting  (PEC)  surface  of  an 
antenna  or  scatterer  is  zero.  Then  the  total  field  E  is  the 
sum  of  the  incident  and  scattered  fields.  This  can  be 
expressed  as 


■^tan  =  ^tan  + 


0  on  S 


(l) 


where  S  is  the  conducting  surface  of  the  antenna  or  scatterer. 
The  subscript  'tan'  indicates  tangential  components.  The 
incident  field  that  impinges  on  the  surface  S  induces  an 
electric  current  density  Js,  which  in  turn  radiates  the 
scattered  field. 

The  scattered  field  can  be  expressed  in  terms  of  a  vector 
potential  A  and  a  scalar  potential  t  as 


where 


E*  (r)  =  -j u>A- V<£ 


(2) 


e  ~jkR 

4  nR 


ds' 


(3) 
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r 


R  =  |r-r'|  (4> 

where  r1  is  the  point  source  location,  and  r  is  the 
observation  point  location  and 


eJJs  M  in R 


(5) 


The  continuity  equation  relates  the  current  and  charge  on  the 
surface 


or 


V-Js  =-jup. 


_2_v-  A  . 

Ci e 


(6) 


(7) 


Substitution  of  (3)  and  (7)  into  (2)  gives  the  following 
equation 


E"  (r)  =  Jon  f  (  JB  (r/)  e  J  ds'  ■*  -2-  V IV- ff  JB  (. r ')  e  J ds']  . 

JJs  4  7t  R  we  JJ  e  inR 

(8) 

This  is  the  most  general  form  of  the  EFIE.  Let  the  free  space 
Green's  function  be  denoted  as 

G(r.r')  -SSM-  a'jk)z'z'}  <’) 

4nR  4n  \  r-r'\ 

so  that  (8)  becomes 


Bs  (r)  =  ju\ij  j  JB(r')  G(r,  r')  ds'  +  -2-V[V-ff  Jt(r')  G(r,  r')  ds'  . 

(10) 

Js  only  depends  on  the  primed  quantities  and  G  is  a  scalar 
that  depends  on  both  the  primed  and  unprimed  quantities.  J,  is 
the  unknown  quantity  to  be  solved  for. 
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For  convinience,  define  an  operator  $£  which  operates  on  Js 

{ [  Jt(r')  G(r,  r')  ds' + -i-V  [V- [  [  J,{r')  G{z ,  r')  ds']  . 

J  J  s  G)€  J  J  s 

<ll> 

After  using  vector  identities  and  modifying  the  last  term 

S£(*7,)  =  jonff  JM(r')  G(r,r')  ds'  +  V[ff  VG(r,r')  ds'  .(12> 

Since  Vg  =  -  V'G,  where  V'  is  defined  with  respect  to  the 
source  quantities, 

ff  VG(r,  z')  •JM(z')  ds'  =  -fj  Vg(z,z')  ‘J.(r')  ds'  .  (13) 

Using  vector  identities  and  the  surface  divergence  theorem,  it 
can  be  shown  that 

2(J.)  =  [(  G(r,r')  - -2- V  ( (V-J,(r')  )  G(r,r/)]]ds/  . 

J  J  s  0)€ 

(14) 

The  EFIE  can  now  be  written  consisely  as 

S£(J,)  =JFtia  .  (15) 

More  detail  on  the  development  of  the  EFIE  is  available  in 
[Ref.  7]  and  in  [Ref.  8]. 

B.  SOLUTION  OF  THE  EFIE  USING  THE  METHOD  OF  MOMENTS 

To  solve  (15)  the  method  of  moments  is  used.  The  method  of 
moments  generates  a  generalized  impedance  matrix  for  the  body 
of  interest.  The  excitation  of  the  body  is  represented  by  a 
voltage  matrix,  and  the  resultant  current  on  the  body  is 
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represented  by  a  current  matrix  [Ref.  2],  Equation  (15)  is 
reduced  to  a  matrix  equation,  which  can  then  be  solved  using 
standard  methods. 

The  analysis  and  solution  presented  here  are  for  the 
special  case  of  bodies  which  are  rotationally  symmetric  about 
an  axis.  Because  of  this  symmetry  a  Fourier  series  expansion 
in  the  angle  of  rotational  symmetry  reduces  the  problem  to  a 
system  of  independent  modes.  This  is  important  from  the  view 
of  computation;  several  small  matrices  must  be  inverted 
instead  of  one  large  one. 

The  first  step  in  the  MM  procedure  is  to  expand  Js  into  a 
series  N 

(16) 

j-i 

where  1^  are  the  expansion  coefficients  to  be  determined  and 
the  are  called  basis  or  expansion  functions. 

The  Jj  are  usually  chosen  to  be: 

(1)  consistent  with  the  behavior  of  the  actual  currat 
on  the  surface, 

(2)  mathematically  convenient  (i.e.  have  certain 
desirable  mathematical  properties  which  will  be 
discussed  later) . 

Equation  (16)  is  substituted  into  (15) ,  which,  because  of  the 
linearity  of  Sf,  becomes 

£  =  Et{J1  (17) 

j 
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Similarity  a  set  of  N  testing  function  {WJ  is  defined. 


Although 

they  need 

not  be  related 

to 

the 

set  of  basis 

functions 

,  they  are 

usually  chosen 

to 

be 

their  complex 

conjugate 

(Galerkin ' 

s  Method) .  Now 

the 

inner  product  is 

defined 

=  W  Jds 

(18) 

where  W  and  J  are  tangential  vectors  on  S.  Taking  the  inner 
product  of  both  sides  of  the  EFIE  (equation  (17)) 

£  (Jj)  )  =  {wi,Ei)  i  =  l,2,3,...,N  .  (19) 

j 

The  subscript  'tan'  has  been  dropped  from  E1  because  the  inner 
product  involves  only  tangential  components.  Now  define 

V^iWi.E1)  (20) 

Zjj  =  (W  j ,  ( J  j ) )  (21) 

so  that  (19)  becomes 

Y/ZijIrVi  i  =  l  ,2  ,  .  .  .  ,  N  .  (22) 

j 

Since  there  are  N  of  these  equations  they  can  be  written  in 
matrix  form 

[Z]  [J]  =  [Vi  (23) 

[Z]  is  the  generalized  impedance  matrix  and  [V]  is  the 
excitation  vector.  The  vector  of  unknown  current  coefficients 
can  be  solved  for  by  inverting  [Z] 


[I]  =  [Z]  1  [ V]  =  [Y]  [V]  . 


(24) 


The  impedance  elements  are 

Zu  =  JJw-  ( j  wAj  +  Vfj)  ds  .  {25) 

S 

If  W  is  thought  of  as  a  current,  the  charge  associated  with  it 
is 


Pi  = 


Now  (25)  can  be  written  as 


(26) 


=  ds  (27) 

S 

or  jkR 

zi}-S!ds'!!ds[  ub) 

S  S  j 

where  primed  donates  a  source  point  (associated  with  Jj)  and 
unprimed  a  field  point  (associated  with  W- ) 


C.  MM  SOLUTION  FOR  THE  SPECIAL  CASE  OF  BODIES  OF  REVOLUTION 

So  far  the  discussion  has  been  for  an  arbitrary  conducting 
body.  In  this  section  we  restrict  the  surfaces  to  those 
generated  by  revolving  a  plane  curve  about  the  z  axis.  Hence 
an  appropriate  set  of  basis  functions  are 


Jmj  =  ucfjt(t)  eJm* 
4j  =  ucff(C)  ejn+ 


m=  0  ,  ±1 ,  ±2  ,  .  .  .  ,  oo 


m=  0  ,  ±1 ,  ±2  , 


(29) 

(30) 
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where  ut  and  are  t-directed  and  (^-directed  unit  vectors.  An 
example  of  body  of  revolution  and  the  coordinate  system  is 
shown  in  Figure  4.  Following  the  method  of  moments,  the 
electric  current  J  on  S  is  given  by 

J-=£  +  .  (31) 

m.j 

For  testing  functions,  Galerkin's  method  is  used  (Wk  =  Jk*)  so 

W*.i  =  utfi(t)e-jl*  (32) 

e-jn*  (33) 

The  Wn  are  orthogonal  to  Jm,  m  *  n,  over  0  to  2n  on  <p ,  and  also 
to  st1  ( JJ  .  Hence,  all  impedance  elements  are  zero  except  those 
for  which  m  =  n,  and  each  mode  can  be  treated  separately.  This 
is  the  major  simplification  introduced  by  taking  advantage  of 
the  rotational  symmetry  of  the  body.  Although  the  sum  over  m 
in  equation  (31) ,  has  an  infinite  number  of  terms,  it  coverges 
rapidly  in  most  cases  and  there-fore  can  be  truncated  at  some 
value  M. 


Figure  4:  Body  of  revolution 
and  coordinate  system 
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The  use  of  (29)  ,  (30)  ,  (32)  and  (33)  to  evaluate  the 


elements  of  (27)  results  in  the  matrix  equation 


'[z“] 

IZ*}' 

t In  ) 

ivf)' 

[ztc] 

[Z**] 

.  ^  . 

Here  the  elements  of  the  submatrices  of  Z  are 


(Z“) 

(z?)  ..  =  (wZ1,'£(J*j)) 

(z*‘  = 

(z^)iJ.  =  (w*i>S£(  J*j)) 

the  elements  of  the  submatrices  V  are 


(36) 

<  V*)  ,  =(«*„.,  B‘> 

and  the  elements  of  the  submatrices  of  I  are  the  coefficients 
in  expansion  (31).  The  solution  to  (34)  can  be  written  as 


'u„£] 

[ 4i 

The  submatrices  of  Y  are  in  general  obtained  after  inversion 
of  the  entire  Z  matrix,  and  not  the  inverses  of  the 
corresponding  Z  submatrices.  The  -n  mode  matrices  are  related 
to  the  +n  mode  matrices  by 


’[ Yntc ]  [Y*]' 

[Yf]  C  Y*] 

[v*] 

[Y“] 

lY-Hi 

’  [y“]  [-y*]' 

[Y^] 

[Y*] 

[-4f]  [Y*] 

Hence,  only  the  n  >  0  mode  matrices  need  be  inverted. 
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D.  MEASUREMENT  MATRICES 

Obtaining  the  scattered  field  from  the  current  J  on  the 
body  S  can  be  done  using  a  measurement  vector 


measurement 


■  //*'■• 


where  Er  is  the  plane  wave  field  scattered  by  a  differential 
current  element  (a  known  function) ,  which  is  weighted  by  the 
current  at  every  point  on  S.  Since  the  current  is  given  by 
(31),  (39)  reduces  to 


measurement  =  [fi]  [J] 


where  R  is  a  measurement  row  vector 


[R]  =  [(ujtE*)]  (41) 

and  I  is  the  current  vector.  Substitution  of  (24)  into  (40) 


gives 


measurement  =  [tf]  [F]  [V]  . 


For  bodies  of  revolution,  the  expansion  for  J  can  be 
separated  into  t-  and  0-directed  components.  Equation  (41) 


becomes 


<Oi  =  Wi'E1) 

(Rn)i  = 


In  partitioned  form, 


,  *  ft*"] 

measurement  =  [[)?„]  [/?„]  ..  . 

[Ff]  [V*] 
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A  special  case  is  that  of  the  radiation  field  measurement. 
Now  the  u-polarized  scattered  far-field  from  the  current  J  on 


S  is 

Ell  =  e-jks  [j?]  [j]  (45) 

AnR 

where  the  elements  of  [R]  are  given  by  (43)  with 

E1  =  uejkr  .  (46) 

This  is  a  unit  plane  wave  with  polarization  vector  u  and 
propagation  vector  k.  An  arbitrary  plane  wave  is  a 
superposition  of  two  orthogonal  components,  for  example  Ee  and 
E  .  Hence,  (46)  can  be  treated  as  two  cases,  one  for  u  =  ue  and 
the  other  for  u  =  u^.  To  distinguish  between  the  two  cases  let 


<Oi  =  (Jni,E0r) 
(i?f)i  =  (J *irEj) 


for  the  0-polarized 


case,  and 
(R?)i  = 

(R?)i  =  Uli.Ef) 


(47) 


(48) 


for  the  0-polarized  case. 
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III.  METHOD  OF  MOMENTS  SOLUTION  FOR  THE  ADE  ANTENNA 


A.  CALCULATION  OF  IMPEDANCE  ELEMENTS 


The  mathematical  analysis  in  this  section  is  based  on 
[Ref.  9:p  14-22]  and  [Ref.  10:p  3-18].  The  expansion  functions 
Jnj‘  and  J  in  equation  (32),  are  chosen  to  be 


u 


Tj(t) 

P 


e  -jn* 


j  =  1,2,  ...  .  NP-2 
n=  0  ,  ±1,  ±2,  .  .  ±<» 


(49) 


U * 


pj  ( t) 

P  i 


ejn+ 


j  =  1,2,...,  NP-1 
n  =  0  ,  ±1  ,  ±2  ,  .  .  .  ,  ±«> 


(50) 


where:  ut  ,  -  unit  vectors  in  the  t  and  <p  direction 

n  -  the  mode  number 

p  -  the  distance  from  the  axis  of  symmetry  to  a  point 
on  the  body  of  revolution 

p-  -  the  value  of  p  at  t=t^  where  t-  is  the  center 
point  of  the  domain  of  the  pulse 

Tj(t)  -  triangle  function  shown  in  Figure  5,  which 
extends  over  segments  j  and  j+1 

Pj(t)  -  pulse  function  shown  in  Figure  6,  which 
extends  only  over  a  single  segment  j . 

The  purpose  of  the  scale  factor  1/p.  is  to  give  (50)  the  same 
dimension  as  (49).  The  derivative  of  Tj(t),  which  is  shown  in 
Figure  8,  is  also  required.  The  mathematical  descriptions  of 
the  functions  are 
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In  Figures  5,6  and  7,  t  is  the  arc  length  along  the  body  of 
revolution  generating  curve  and  the  primes  denote 
differentiation  with  respect  to  the  argument.  It  is  assumed 
that  the  generating  curve  consists  of  NP-1  straight  line 
segments  where  NP  is  an  integer  greater  than  or  equal  to  3. 
The  j -th  segment  extends  from  t.'  to  t  7+1  and  has  length 


For  testing  functions  Galerkin's  method  is  used 


"ai  =  iJni) 


(55) 


(56) 


where  *  signifies  the  complex  conjugate. 


In  terms  of  the  basis  functions,  the  current  is  expressed 


»  NP-2  NP- 1 

n*-»  i*l  j-l 


After  applying  the  method  of  moments  testing  procedure  a 
matrix  equation  of  the  form  of  equation  (25)  results 


I  =  Z''V 


where 


I  is  the  vector  of  current  coefficients,  V  is  the  excitation 
vector  and  Z  is  the  method  of  moments  impedance  matrix. 

The  impedance  matrix  can  be  partitioned  according  to  the 
current  and  testing  function  components.  For  each  azimuthal 


mode 


ytt  yt* 

^ n  6 n 
^ n  & n 


For  convenience  let 
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MT  =  NP-2  (  number  of  surface  triangles  )  .... 

\  ®  ^  ) 

MF  =  NP- 1  (  number  of  surface  pulses  ) 

Thus  the  dimensions  of  the  submatrices  of  Z  are 

Z“  -  MTxMT 
Z„*  -  MTxMP 

(62) 

-  MPxMT 
Z **  -  MPxMP 


Althought  n  can  take  on  values  between  ±°o,  the  sum  must  be 
truncated  to  some  finite  value,  M.  When  multiple  modes  are 
used,  the  complete  impedance  matrix  becomes 


Fortunately  this  a  block  diagonal  matrix;  the  inverse  of  Z  is 
also  a  block  diagonal  matrix  with  diagonal  submatrices  (Z  ss)'1. 

To  explicitly  evaluate  the  impedance  elements,  the  basis 
and  testing  functions  are  used  in  equation  (28) 


(64) 
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where  p  and  q  signify  all  possible  combinations  of  t  and  <p , 
and  T]  is  the  intrinsic  impedance.  It  is  only  necessary  to 
evaluate  a  few  of  these;  the  remainder  can  be  obtained  by  the 
following  mode  symmetry 

[Zfl“]  ij  [z.%  _ 

[4%  [Z*  }ij 

and  symmetry  between  positive  , 

[Z*]  2j 

[z'tiij  [Z*]^  "  - 
For  a  body  of  revolution 

kp 

//ds  =  / p( 

s  0 

and  similarly  for  the  primed  integration.  The  dot  products  are 

ff(  t)  f/(  t')  Up‘llq  .  (68) 

The  functions  f;  and  fj  are  either  the  pulse  or  triangle 
depending  on  p  or  q. 

The  unit  vector  dot  products  can  be  defined  in  terms  of 
integration  variables  in  Figure  4, 

Uj-  ut  =  sin  (v)  sin  (v7 )  cos  ( <p-<p ' )  +  cos  (v)  cos  (v7 ) 

ut'u#  =  ~sin(v/)  sin(0-07)  (69) 

u^-  ut  =  sin  (v)  sin  (0-07 ) 
n'$-  u,  =  cos  (<p-<p') 


[Z^lji  -[4% 
-[Z?]u 

and  negative  modes 


(65) 


[4%  [Z^ij 


(66) 


2* 


<J>)  dt  J d<$> 


(67) 
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Here  v  is  the  angle  between  the  t  direction  and  the  Z-axis, 
being  positive  if  ut  points  away  from  the  Z-axis  and  negative 
if  ut  points  toward  it.  Using  equations  (51),  (52),  (53)  and 
(54) ,  the  elements  of  Z  are  given  by 


(z“)  lj  =  j  j  dt  j  dt'  (kTi{t)  Tj  { t')  (G1sin  ( v)  sin  ( v') 

tj  tj 

+  G3  cos  ( v)  cos  ( v7)  )  -  G3  T  ( t)  Tj  ( t')  ) 

3  3  dt  1  dt'  3 


(70) 


(4£> 


p. 


1  ^j*l 

|  dtPjit)  f  dt'  (kpTj(t')  G2sin(v') 


(71) 


*nGy—.TAt') 
3  dt'  J 


ci«-2  ^;*l 

(£„*)„  =  y:  /dt/  dt'Pjit')  ikp'Tjit)  Gasin(v) 


t/  tj 


(72) 


tj.l 


= 


J 


P.PJ 


f  dt  Pjit)  f  dt'Pj(t')  (kpp'  G,  -  nz  G3) 


(73) 


where 


G0  =  2  j”  d<J> 


e~jkR  ,  4> 


,  sin2  ( )  cos  (n<{>) 
icl?  2 


(74) 


/p  -  j*/? 

d<{>  — - —  cos  (<|>)  cos  (n<f>) 
kR 


(75) 
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G2  =  fd<b  sin  (4>)  sin  (jq4>)  (76) 

0 

* 

r  a  ~jkR 

G3  =  G0  +  G1  =  fd<t>  e .  — COS  (r?4>)  (77) 

i?  = 

The  parameters  p,  z,  and  v  depend  on  t,  while  p',  z '  , and  v' 
depend  on  t '  . 

To  reduce  the  number  of  numerical  integrations,  the  6- 
function  approximation  is  used  for  the  integral  with  respect 
to  t  in  (70)  -  (73) 

tp.i 

I  f{t)  dt  -  f(tp)  Ap  (79) 

t ; 

where  f(t)  is  the  relevant  integrand  and,  as  indicated  in 
Figures  6  and  7, 


tp 

~  u p+  tptl) 

(80) 

*p  - 

(  1  —  )  • 

(81) 

Application  of  (79)  to  (70)  -  (73),  leaves  only  one 

integral  with  respect  to  t  in  each  equation.  These  are 
evaluated  numerically  as  are  the  integrals  with  respect  to  <p 
in  the  G  functions  defined  earlier.  Another  integration 
shortcut  is  achieved  by  noting  that  when  testing  is  performed, 
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the  integrals  in  adjacent  impedance  elements  that  span  a 
common  segment  are  the  opposite  of  each  other.  Thus,  for  the 
impedance  elements 


V'  - 

(££)  a  =  J*pj  dt'  (kti{tp)  Tjit')  (Gxsin  ( vp)  sin(  v') 

^ Q 

+  G3  cos  (  v  )  cos  (  v1)  )  -G3  [  Ti  ( t)  ]  t  — —  Tj  ( t') 

3  P  3  dt  2  Cp  dt,  J 


(82) 


(zT)PJ  - 


-A_P_(tr)  f  dt'  (kTAt1)  G2  sin(v')  +  -£L  G3  — —  T  •  ( t’) 
p  p  p  J  J  2  PP  3  dt'  J 


(83) 


2?)iq  =  *p/  dt'PQ(t')  (Ae_ri(£p)G2sin(vp)+-^G3[^ri(t)]tp) 

r  P  M  -nr 


(84) 


^<7*1 


a,)  . 

V  q  P pP  gr 


(85) 


The  A  denotes  the  contribution  due  to  integration  over  the 
region  A^,  and  therefore  i  is  p-1  or  p  and  j  is  q-1  or  q.  Note 
that  all  the  unprimed  quantities  that  depend  on  t  are  now 
evaluated  at  t=t  due  to  the  approximate  integration  method. 
Further  substitution  for  the  functions  P(t)  and  T(t)  give  the 
final  equations 
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(86) 


ic  A„  A. 


(zZc)ij  =  j — I-2-)  <Glasin(  vp)  sin(vq)  +G3acos  (  vp)  c  os(v,)) 
( Gljbsin  (  vp)  sin  (vQ)  +G3bcos(vp)  cos  (vg) 


u  -  8 


8 

- o - JG2t 


+  (^'  G3  ) 

2  PD  3a 


*V^(y  )  G2>_  )  G;e 


(87) 


<**>„-  (^4^)  +  (.1}p-i  (^) 


2  Pc 


(88) 


(Opq  =  2j[( 


<G„  + 


Aqsin(vq) 

2pl 


nA„  nA n 

v-WW*  (89> 


Gj.  = 


C<7*  1 

J  Gm  ( t'  -  tq)  dt'  w  =  0,1,2 


(90) 


<7*1 


7J£> 


(  —  )2  f  (t'-  tq)  Gm(t7  -  tq)  dt1  m  =  0,1,2 

A  „  J 


(91) 


Substitution  of  Rp  for  R  in  (75)  -  (77)  gives 


G0(  t'- tq)  =  2  fd<t> 


- jkRp  x 

sin2  ( )  cos  (n<J>) 


c  kRP 


(92) 


~jkF 


Gl{t'-tq)  =  J  d(j>  -ekR  P  cos  (<j>)  cos  (n<{>) 


kR 

0  P 


(93) 
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(94) 


Evaluation  of  the  integrals  in  (92)  through  (94  )  is  done 
using  nt~point  Gaussian  quadrature 

b  nt 

ff(x)  dx  -  c.f(aj  <«> 

a 


where  cm  are  the  weights  and  xm  are  the  zeros  of  the  Legendre 
polynomials.  Typically  nt  =  2  is  sufficient  for  the  t'- 

integration.  For  the  0-integration  the  appropriate  value 
depends  on  the  maximum  radial  extent  of  the  surface.  (An 
approximate  rule  of  thumb  is  10  points  per  wavelength  of 
circumference.)  The  calculated  values  of  Gm  are  substituted 
into  (92)  and  (93)  in  order  to  evaluate  Gm  and  G,.  .  The 
resulting  values  of  G^  and  G^  are  then  substituted  into  (86) 
-  (89)  to  evaluate  the  elements  of  the  impedance  matrix. 

The  general  MM  solution  can  be  applied  to  multiple 
detached  surfaces  as  in  the  ADE  antenna.  The  extension  to 
multiple  surfaces  is  primarily  a  bookkeeping  problem  that 
complicates  the  computer  code.  A  simple  way  of  calculating  the 
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current  on  multiple  bodies  is  to  perform  the  calculation  as  if 
they  were  attached  and  then  set  the  appropriate  impedance  and 
excitation  elements  to  zero. 

Figure  8  shows  the  two  reflector  surfaces  of  the  ADE 
antenna  connected  by  a  fictitious  "dummy  segment."  After  the 
Z-matrix  calculation,  any  element  that  spans  the  mth  segment 
is  set  to  zero,  as  shown  in  Figure  9  for  the  case  of 
triangular  basis  functions  and  Figure  10  for  pulses. 


1*  1  2  3  n-1  b  *5  If- 2 


Figure  9:  Triangle  segments  Tj(t) 

To  avoid  a  singular  matrix  the  diagonals  of  the  zeroed  rows 
and  columns  are  set  to  1.  The  final  matrix  structure  is 
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illustrared  in  Figure  11.  This  method  does  not  use  computer 
memory  efficiently,  but  is  suitable  for  antennas  with  only  a 
few  detached  surfaces. 


Figure  10:  Pulse  segments  -  P.  (t) 


2 

\ 

m 

\J 

v‘ 

1 

m-1 

z* 

2 

]♦ 

v* 

2m 

1  on  the  diagonal 
nulled  elements 


Figure  11:  Resulting  impedance  matrix  for  multiple 
scatterers 

B.  EVALUATION  OF  THE  EXCITATION  VECTOR 

The  excitation  vector  for  a  0  or  ^-polarized  wave  incident 
in  the  cp- 0  plane  is  defined  by 


krje 


-jkt-r 


(98) 


where 


kt  =  -  k  (uxsin  (0t)  +  u,sin  (0t)  ) 


(99) 
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(100) 


u8  =  uxcos(0t)  -uzsin(0t)  .  (101) 

0t  is  the  angle  of  incidence  and  ux,  uy>  and  uz  are  unit 
vectors  in  the  x,  y,  and  z  directions,  r  is  the  position 
vector  from  the  origin  to  the  test  point  on  the  body.  The 
elements  Vnr'*  are  given  by 

V&*  =  j-ffw^-Eids  .  (102) 

A  detailed  development  of  the  excitation  vector  can  be  found 
in  [Ref.  10:  p  27-31].  The  final  results  are  obtained  by 
substituting  (98)  -  (101)  in  (102) .  Using  the  same  notation 

developed  earlier 


to  jn‘VkA  sin(v  )  cos  (0t) 


j"7TkApcos(vp)sin(0t) 

vni  “  - 7 -  5 - *  r 


_  ; ,  jn+17rkADsin  (v_)  cos  (0.) 

+  5 - —  (Fn.1tb-Fn.ltt 


jn7rkA  cos(v  )sin(0t) 


'  nb 


(103) 


j‘*k*pc 08  (6t)  .  ipsin(Vp) 

vnp  0  [  [rn*l,a  +  rn-l.aJ  o  [rn*l,b  rn-l.b>> 


2  P, 


(104) 


-  in'lnkAp 
**  2 


(iF  -  F  )  .  APSin(vP}  ,F  _F  ))<105> 
urn*l,a  cn-l.a>  ?  ('rn*l.b  rn-l,b>> 

*  I3  D 
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(106) 


j  nnkApsin(vF)  ^  ^  N  (-1)P'JJ  /cApsin  ( vp) 

^  n+ 1 ,  a  +  * n-1.  a  '  ^ 

•(Fn.l,b+Fn-i.b>  • 


New  functions  have  been  defined  in  (103)  -  (106), 


Fma  =  —  f  Jm(kp  sin  (<t>c)  )  e-J**co9(+d  ;  m=n-l,n,n+l  ^107^ 

A  w 


tp-i 


Fmb  '  (  —  ) 2  /*  (  £  -  t)  J*  Up  sin  (0,)  )  eJ 

A  _  J 


kz coe(6,!  ..  _  „  -  _  „  , 

-  at  ;  m-  n-1 ,  n,  ml 


(108) 


Jm (kpsin  (0t)  )  is  the  m-th  order  Bessel  function,  which  results 
from  the  integration  in  0,  and 


p  =  pp+  ( t- tp)  sin(vp)  (109) 

z  =  zp  +  ( t-  tp)  cos  ( vp)  .  (110) 

The  remaining  integrations  are  evaluated  using  a  Gaussian 
quadrature  formula  as  in  the  case  of  the  impedance  elements. 

At  this  point  all  of  the  quantities  required  to  solve  the 
method  of  moments  matrix  equation  have  been  defined.  An 
existing  computer  code,  originally  developed  by  Mautz  and 
Harrington,  was  modified  to  compute  these  elements  for  the  ADE 
antenna.  The  computer  code  operation  is  described  in  the 
following  chapter. 
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IV.  COMPUTER  PROGRAM 

In  this  section,  the  FORTRAN  program  which  implements  the 
numerical  solution  of  CHAPTER  III  is  discussed  (the  computer 
listing  is  given  in  Appendix  A)  .  The  program  consists  of  a 
MAIN  PROGRAM  and  the  subroutines  ZMAT ,  SPHERE,  PLANE,  PLAN, 
DECOMP,  SOLVE  and  the  functions  BLOG,  ETHETA,  EPHI . 

The  main  program  obtains  the  radiation  pattern  and  gain  of 
an  ADE  antenna  for  an  ideal  feed  at  body  of  revolution 
coordinate  system  origin.  The  subroutine  SPHERE  calculates  the 
elements  of  excitation  vector  in  (23)  for  spherical  wave 
incidence.  The  subroutine  ZMAT  calculates  the  elements  of  the 
moment  matrix  in  (24).  The  function  BLOG  is  called  by  ZMAT. 
The  subroutines  DECOMP  and  SOLVE  solve  the  matrix  equation 
(24)  for  Int,#.  The  subroutine  PLANE  calculates  the  plane  wave 
receive  vector  elements  need  to  obtain  the  radiation  field. 

A.  THE  MAIN  PROGRAM 

The  main  program  is  divided  to  two  main  parts.  The  ADE 
ANTENNA  RADIATION  PATTERN  (lines  1  -  334)  and  the  GAIN 
PROGRAM ( lines  335  -  453).  Both  parts  are  based  on  Mautz  and 
Harrington’s  computer  program  for  bodies  of  revolution. 

1.  The  ADE  Antenna  Radiation  Pattern 

This  portion  of  the  main  program  computes  the 
radiation  pattern  of  an  ADE  type  of  antenna.  An  ideal  feed  (a 
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point  source  with  a  directive  pattern)  is  at  BOR  coordinate 
system  origin  and  therefore  only  n  =  ±1  are  required.  Lines  17 
73  are  basic  initializations  including  integration 
coefficients  for  the  ^-integration  (NPHI)  ir.  lines  31  -  37, 
and  integration  points  for  the  t-integration  (NT)  in  lines  43 
-  48.  The  antenna  geometry  is  defined  in  lines  49  -  59:  the 
total  number  of  generating  points  (NP) ,  the  number  of  the  main 
reflector  points  (NP1)  and  the  number  of  the  subreflector 
points  (NP2) .  Lines  60  -  73  define  the  feed  pattern  which  is 
of  the  form 


~  7  k  R  4*  a 

Et»ad  =  — p —  (Qcosn(8)  cos(4>)  +  ^cos'MB)  sin(4>))  (111) 

R 

where  n  =  FEXP  and  m  =  FHXP  in  the  code.  The  simple  case  of  a 
Lambertian  source  is  n,m  =  1. 

The  next  portion  (74  -  211)  of  this  part  cf  the  program 
generates  an  array  of  surface  contour  pc ' nts  for  calculating 
the  moment  matrix  and  the  excitation  vector.  Lines  74  -  113 
define  the  antenna  parameters  which  are  discussed  in  Chapter 
I.  To  create  the  contour  of  the  main  reflector  (which  is 
paraboloid)  ,  1J1  is  varied  in  the  equation 


R. 


2  f 


(112) 


1  +  cos  (\jl ) 

t|i  is  the  polar  angle  from  the  paraboloid  vertex  and  Rm  is  the 
distamce  from  the  focus  to  the  surface  point  at  angle  i|».  The 
coordinates  of  the  surface  are 
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p  =  R^cos  ( tj/) 
z  =  R^sin  (ijj)  . 


(113) 


For  the  subreflector  (which  is  ellipsoid) , 

p  =  A4  ( 1  -  ECC* ) _  .  (114) 

s  1  ±  ECC-c os(0) 

The  values  of  kp  and  kz  are  stored  in  the  arrays  RH  and  ZH. 


After  ZH  and  RH  are  determined,  the  program  calls  the  SPHERE 
subroutine  to  calculate  the  excitation  coefficients  and  ZMAT 
to  get  the  moment  matrix.  Next,  lines  212  -  226  take  care  of 
the  "dummy  segment"  between  multiple  surfaces,  mentioned  in 
th®  previous  chapter.  In  1  es  227-  240  the  subroutines  DECOMP 
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and 


and  SOLVE  solve  the  ma+  rix  equation  (24)  for  Inl  and  In*, 
in  the  lines  241  -  334  the  far  field  patterns  in  the  E  and  H 
planes  are  computed  according  to  equations  (40)  -  (49). 

2.  Gain  Program 

The  second  half  of  the  main  program  calculates  the 
gain  and  the  efficiency  of  the  ADE  antenna.  Far  field  pattern 
integration,  based  on  the  previously  calculated  currents,  are 
used  to  get  the  total  radiated  power.  Lines  3  45  -  360  are 
initialization  and  specifying  integration  coefficients  for  the 
0  (NPHI )  and  t  (NT)  integrations.  A  normalization  is  loop 
incorporated  in  lines  351  -  382  to  obtain  E^.  In  lines  383  - 
402  the  0  and  0-integration  intervals  are  defined 
The  directivity  is  computed  from  the  formula 

2  It  Jl 

\2R2SIN(Q)  d0a'(j)  (115) 

0  0 


DIRECTIVITY  = 


4*  [iwl2 

0* 


(116) 


where  E(0,0)  is  the  total  electric  field  (i.e.,  that  radiated 
by  the  feed  plus  that  scattered  from  the  reflector)  .  Emax  is 
obtained  from  the  normalization  loop. 


B.  THE  SUBROUTINE  SPHERE 

The  subroutine  SPHERE  computes  the  spherical  wave 
excitation  vector  for  the  body  of  revolution.  The  source  is 
assumed  to  be  located  at  the  origin  of  the  BOR  coordinate 
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system.  The  arguments  are  defined  in  the  program.  It  returns 
R,  a  column  vector  of  length  N  containing  V*  and  V*  for  mode 
n  =  +1,  in  that  order. 

C.  THE  SUBROUTINE  2 MAT 

The  subroutine  ZMAT  calculates  the  moment  matrix  in  (23) 
for  n  =  1  and  stores  them  in  Z,  which  is  the  only  output 
argument.  The  storage  of  the  Zn  submatrices  in  Z  is  as 
follows: 

(Z„“)  e  z(i  +  N*  (j-i) ) 

(zf )ij  6  Z(i+N*  (j-1)  *  NP-2) 

(Z^)  ij  €  Z  (i  +  N*  (j-1)  +  (NP-2)  *N) 

(Z**)  i:j  €  Z(i+N*  (j-1)  *  (NP-2)  *N+ NP-2) 

where  N=2*NP-3.  More  details  about  the  ZMAT  subroutine  can  be 
found  in  [Ref.  10:  p.  43-56]. 

D.  THE  SUBROUTINES  DECOMP  AND  SOLVE 

The  subroutines  DECOMP  and  SOLVE  solve  a  system  of  N 
linear  equation  in  N  unknowns.  The  input  to  DECOMP  consist  of 
N  and  the  N  by  N  matrix  of  coefficients  on  the  left-hand  side 
of  the  matrix  equation  stored  by  columns  in  Z.  The  output  from 
DECOMP  is  IPS  and  Z.  The  output  is  fed  into  SOLVE.  The  rest  of 
the  input  to  SOLVE  consists  of  N  and  the  column  of 
coefficients  on  the  right-hand  side  of  the  matrix  equation 
stored  in  B.  The  solution  to  the  matrix  equation  is  in  X. 
Further  details  concerning  DECOMP  and  SOLVE  can  be  found  in 
[Ref.  13:  p.  46  -  49], 
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E.  THE  SUBROUTINES  PLANE  AND  PLAN 

The  subroutines  PLAN  and  PLANE  are  essentially  the  same 
and  compute  the  plane  wave  receive  and  excitation  vectors  for 
the  BOR  according  to  (103)  -  (106).  They  are  stored  in  R  which 
is  the  only  output  argument.  The  rest  of  the  arguments  of 
PLANE  and  PLAN  are  input  arguments  and  have  all  been  defined 
previously,  except  NF  which  is  always  1  as  implemented  here. 
The  receive  vectors  are  returned  in  the  following  order:  Rte, 
R*6,  Rttf,  Rw.  These  subroutines  were  written  by  Mautz  and 
Harrington,  and  are  described  in  detail  in  [Ref.  12:  p.  57  - 
62]  . 
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V 


COMPUTED  RESULTS 


An  evaluation  of  practical  ADE  antenna  configurations  has  been 
performed  using  the  MM  computer  code.  The  main  reflector 
diameter  was  6"  at  44  GHZ  (22.352  wavelength),  which  was 
chosen  the  same  as  the  M. I . T-LINCOLN  Lab  design  [Ref.  6].  A 
second  fixed  parameter  was  the  focal  length  to  main  reflector 
diameter  ratio  (FOD)  which  was  0.15.  The  displacement  distance 
from  the  axis  of  symmetry  to  inside  edge  of  the  main  reflector 
is  (d) .  The  ratio  of  (d)  to  the  main  reflector  diameter  (dOD) 
was  fixed  at  0.2.  The  edge  angle  of  the  paraboloid  (ALPHA)  is 
equal  to  105 1  The  subreflector  edge  angle  and  the 
eccentricity  (or  equivalently,  the  feed  location,  Rf)  were 
varied.  The  range  of  design  parameters  that  was  investigated 
are  summerized  in  Table  1.  For  each  of  these  antennas  the  gain 
was  calculated.  The  gain  was  referred  to  that  of  a  uniformly 
illuminated  unblocked  aperture  to  arrive  at  an  efficiency. 
Because  MM  includes  all  of  the  interactions  between  the 
surfaces,  all  of  the  loss  mechanisms  (edge  diffractions, 
surface  waves,  multiple  reflections,  etc.)  are  considered. 
Only  the  feed  is  assumed  to  be  ideal. 

Three  typical  reflector  configurations  are  shown  in 
Figures  13,  14,  and  15.  The  principal  plane  patterns  are  given 
in  Figures  16,  17  and  18.  The  optimum  occur^d  for  BETA  =  45 1 
For  smaller  values  of  BETA  the  spillover  increased.  For  larger 
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values  of  BETA  the  subreflector  edge  begins  to  turn  back  on 
itself,  which  causes  edge  diffraction  to  increase  and  some 
partial  blocking  of  the  wave  reflected  from  the  paraboloid. 


TABLE  1 


ADE  INVESTIGATED  PARAMETERS 


PARAMETER 

RANGE  OF  VALUES 

DM 

22 . 352 A. 

FOD 

0.15 

dOD 

0.2 

ALPHA 

105° 

BETA 

30°  to  55° 

Rf 

0.251  to  2.751 

The  subreflector  shape  also  affects  the  radiation  in  the 
rear  hemisp.'  ere  (the  backlobe)  .  The  computer  code  assumes  that 
there  is  a  hole  of  diameter  d  in  the  center  of  the  main 
reflector.  According  to  ray  optices,  there  will  be  no  energy 
passing  through  the  hole.  In  this  case  the  subreflector  is 
rather  small  in  terms  of  wavelengths,  so  its  scattered  field 
is  not  very  directive.  Therefore  some  energy  reflected  from 
the  ellipsoid  will  pass  through  the  hole.  The  effect  is  more 
noticeable  when  the  subreflector  projected  area  is  larger,  as 
in  the  case  of  larger  BETAS. 
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Figure  15:  ADE  ANTENNA  SHAPE  FOR  BETA=55° 
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Figure  16:  E-FIELD  RADIATION  PATTERN  FOR  BETA=45° 


Figure  17:  E-FIELD  RADIATION  PATTERN  FOR  BETA=3 0  ° 


Figure  18:  E-FIELD  RADIATION  PATTERN  FOR  BETA=55 


The  effect  of  moving  the  feed  point  along  the  axis  of 
symmetry  was  investigated.  In  practice  it  is  desireable  to 
have  the  feed  as  close  to  the  main  reflector  as  possible.  This 
reduces  the  transmission  line  length  and  also  minimizees  the 
interaction  between  the  feed  and  subreflector.  However,  the 
subreflector  edge  angle  is  smaller  in  this  case  and  there  will 
be  more  spillover.  Spillover  is  the  primary  contributor  in 
reducing  the  gain  of  small  dual  reflector  antennas.  To  reduce 
spillover  requires  a  larger  feed  (narrower  beamwidth)  which  in 
turn  increases  the  interaction  between  the  feed  and 
subreflector.  The  subreflector  will  be  in  the  near  field  of  a 
large  feed,  so  the  illumination  amplitude  and  phase  will  not 
be  uniform.  The  distance  of  the  feeder  from  the  main  reflector 
was  changed  from  0.25  wavelength  to  2.75  wavelength  in 
increments  of  a  quarter  wavelength.  At  the  same  time,  for  each 
value  of  Rf,  the  subreflector  edge  angle  (BETA)  was  changed 
from  30°  to  55°  The  eccentricity  is  directly  related  to  the 
feeder  location  (there  is  one-to-one  correspondence) .  The  goal 
was  to  optimize  the  efficiency  over  all  of  the  parameter 
ranges  described  above.  The  results  are  summerized  in  Figure 
19.  The  highest  efficiency  is  48.3%  for  a  Lambertian  feed 
( COS (0 ) )  ,  when  it  is  located  at  1.5  wavelengths  from  the  main 
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A  second  investigation  was  performed  to  examine  the  effect 
of  the  feed  directivity.  The  shape  of  the  theoretical  feeder 
was  changed  from  COS(0)  to  COS2(0).  The  geometric  parameters 


Figure  19:  SUBREFLECTOR  SHAPING  VERSUS  EFFICIENCY  FOR 
DIFFERENT  FEEDER  LOCATIONS  (COS (6)) 

were  varied  through  the  same  range  as  in  the  case  of  Figure 
19.  As  expected,  by  increasing  the  feed  directivity  the 
efficiency  has  been  improved.  However,  the  optimum 
configuration  remained  the  same.  The  assumed  ideal  feed 
pattern  can  be  integrated  directly  to  find  the  spillover 
energy.  This  was  done  for  the  optimum  configuration  (48.3% 
efficiency)  and  the  spillover  loss  was  to  found  to  be  3  dB. 
Removing  this  loss  from  the  48.3%  gives  96.6%,  showing  that 
the  feed  spillover  is  indeed  the  major  source  of  energy  loss. 
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Figure  20:  SUBREFLECTOR  SHAPING  VS.  EFFICIENCY  FOR 

DIFFERENT  FEEDER  LOCATIONS  (COS2(6)) 


VI.  SUMMARY  AND  CONCLUSIONS 

In  the  previous  chapters  the  EFIE  was  derived  for 
conducting  surfaces  and  then  specialized  to  multiple  surfaces 
of  revolution  with  the  ADE  shape.  A  computer  program  was 
written  to  solve  for  currents,  far  field  radiation  pattern, 
gain  and  efficiency.  Efficiency  as  a  function  of  subreflector 
shape  (BETA)  and  the  feeder  location,  was  investigated  for 
various  feed  directivities. 

The  computed  results  that  were  obtained  were  predicable. 
Maximum  efficiency  for  an  ideal  source  such  as  Lambertian  was 
48.3%  for  6"  diameter  ADE  antenna  at  44  GHZ.  The  first 
sidelobe  level  was  -12  dB  below  the  main  beam.  The  second 
sidelobe  was  -25  dB  down  at  angle  of  9°  from  the  main  lobe, 
which  is  at  the  same  location  measured  in  [Ref.  6].  The 
results  show  how  the  far-out  sidelobes  drop  off  rapidly,  which 
helps  to  reduce  interference  from  outside  sources  which  are 
not  located  close  to  the  main  beam.  The  back  lobe  radiation 
varies  with  the  subreflector  shape  because  of  the  hole  at  the 
center  of  the  main  reflector.  It  is  in  the  -25  dB  range 
relative  to  the  main  lobe  which  is  typical  for  this  small 
reflector  size.  The  efficiency  of  the  ADE  is  higher  than  that 
of  a  Cassegrain  antenna  of  the  same  size  and  feed 
illumination . 
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The  computed  and  measured  patterns  for  the  configuration 
of  [Ref.  6]  are  for  a  dual  band  feed,  which  was  not  included 
in  the  computed  results.  The  program  neglects  feed 
interactions  with  the  reflector  surfaces.  In  [Ref.  6.]  a 
circular  waveguide  is  located  very  close  to  the  subdish.  Hence 
the  results  are  noticeably  different  for  both  the  gain  and 
efficiency.  However,  it  has  been  demonstrated  that  the  program 
gives  an  optimum  point  for  the  antenna  geometry  that  is 
relatively  insensitive  to  the  feed  pattern.  This  provides  a 
starting  point  for  the  antenna  design  process.  Once  a 
particular  feed  is  selected,  a  development  model  can  be  built 
and  optimized  by  tuning  the  feed  to  minimize  the  feed  and 
subreflector  interactions. 

The  next  logical  step  is  to  include  a  more  accurate  feed 
model  into  the  antenna  code  so  that  the  interactions  between 
it  and  the  other  surfaces  are  accounted  for.  Simple  feeds  can 
be  modeled  directly;  more  complicated  feeds  can  be  simulated 
using  measured  data.  Since  the  suoref lector  often  lies  in  the 
near  field  of  the  feed,  a  spherical  wave  expansion 
representation  of  the  measured  feed  field  should  be  used. 

The  use  of  low-blockage  feed  types  such  as  dipoles  and 
parasitic  elements  should  be  investigated.  These  types  of 
feeders  will  minimize  the  interaction  between  the  feed  and 
subdish  and  can  be  included  by  slightly  changing  the  program. 
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APPENDIX 


COMPUTER  CODES 

1  C  *ADE  ANTENNA  RADIATION  PATTERN* 

2  C 

3  C  RADIATION  PATTERN  OF  AN  ADE  ANTENNA 

4  C  BASED  ON  MAUTZ  AND  HARRINGTON'S  COMPUTER  PROG  FOR  BORS . 

5  C  SPHERICAL  SOURCE  IS  AT  BOR  COORDINATE  SYSTEM  ORIGIN. 

6  C  (THIS  IMPLIES  n=+l,-l  ARE  THE  ONLY  TWO  NONZERO  MODES.) 

7  C  VARIABLE  FEED  EXPONENT;  GAIN  AND  SPILLOVER  LOSS  ARE 

8  C  COMPUTED 

9  C 

10  C  ICALC=0  CALC  CURRENTS  AND  FIELD 

11  C  =1  CALC  FIELDS  ONLY  (READ  CURRENTS) 

12  C  ICWRT-0  WRITE  CURRENTS  ON  DISC  FILE 

13  C  ISEG=0  PRINT  SURFACE  SEGMENT  POINTS 

14  C  I PRINT=0  PRINT  TABLE  OF  FAR  FIELD  POINTS  —  OTHERWISE  ONLY 

15  C  GAIN  AND  SPILLOVER  ARE  PRINTED 

16  C 

17  COMPLEX  EP , ET ,Z(300000) , R (4  000 ) ,B(2000) ,C(2000) ,  U ,  UC 

18  COMPLEX  ETT , EPT , ETST, ETSP, EPST , EPSP , EC , EX 

19  DIMENSION  RH ( 2000) , ZH ( 2000 ) , XT ( 10) , AT ( 10 ) , IPS ( 500 ) 

20  DIMENSION  AG ( 100 ) , A ( 300) , X ( 300 ) , EXP ( 12 00 ) , ANG ( 1000 ) 

21  DIMENSION  ECP ( 1000) , XG ( 100) , ECV ( 1000) , EXV ( 1000 ) 

22  DIMENSION  PHC ( 1000 ) , PHX ( 1000 ) , Q ( 3 0 ) , S ( 3 0 ) 

23  DATA  PI, START, STOP, DT/ 3. 14159,0. ,180. ,1./ 

24  DATA  I CALC/ 0/ , ISEG/0/ , IPRINT/ 0/ , ICWRT/0/ 

25  DATA  NT/2/, XT(1) ,AT(1)/. 5773503, 1./ 

26  C  DATA  NT./4/,XT(l)  ,  XT  (  2  )  ,AT(1)  ,AT(2)/.  33998104, 

27  C  *.8611363115, .6521451548, .3478548451/ 

28  C 

29  C  READ  GAUSSIAN  CONSTANTS 

30  C 

3  1  OPEN ( 1 , FILE= '/home/srvl/vered/thesis/outgaus ' , 

32  *STATUS= ' old ' ) 

33  READ ( 1 , * )  NPHI 

34  DO  3  K=1 , NPHI 

35  READ ( 1 , * )  X ( K) , A ( K) 

36  3  CONTINUE 

37  CLOSE ( 1 ) 

38  RAD=PI/180 . 

39  BK=2 . *PI 

40  U=(0.,1.) 

41  U0=(0.,0.) 

42  UC=-U/4 . /PI 

43  NT2=NT/2 

44  DO  1  K=1 , NT2 
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1  K1=NT-K+1 

2  AT (K1 ) =AT ( K) 

3  XT (Kl) =XT (K) 

4  1  XT ( K) =-XT ( K) 

5  C 

6  C  NP1=NUMBER  OF  MAIN  REFLECTOR  PTS . ;  NP2=NUMBER  OF 

7  C  SUBREFLECTOR  PTS. 

8  C 

9  NP1=9 1 

10  NP2=19 

11  NP=NP1+NP2 

12  MP=NP-1 

13  MT=MP-1 

14  N=MT+MP 

15  write(6,*)  'np,n=',np,n 

1 6  C 

17  C  SPECIFY  FEED  PATTERN  WHICH  IN  OUR  CASE  IS  LAMBERTIAN  WITH 

18  C  N=1  AND  IS  LOCATED  AT  THE  ORIGIN. 

19  C 

20  C  FEED  FUNCTION  EXPONENT: 

21  C  ETHETA^COS (THETA) **FEXP*COS (PHI) 

22  C  EPHI=COS (THETA) **FHXP*SIN ( PHI ) 

23  C 

24  FEXP=1 . 0 

25  FE2=2 . *FEXP+1 . 

26  FHXP=1 . 0 

27  FH2=2.*FHXP+1. 

28  DIR=2 . *FE2*FH2/ ( FEXP+FHXP+ 1 . ) 

29  PRAD=1./FE2+1./FH2 

30  C 

31  C  READ  THE  ADE  INPUT  PARAMETERS: 

32  C 

33  C  DM=MAIN  REFLECTOR  DIAMETER  (PARABOLOID) 

34  C  Rf eed=THE  DISTANCE  BETWEEN  THE  FEEDER  AND  THE  MAIN  REFLECTOR 

35  C 

3  6  C  FOD=THE  RATIO  BETWEEN  THE  FOCUS  AND  THE  DIAMETER  OF  THE  MAIN 

37  C  REFLECTOR 

38  C  d/2=DISTANCE  BETWEEN  THE  EDGE  OF  THE  MAIN  REFLECTOR  AND  THE 

39  C  AXIS  OF  SYMMETRY 

40  C  dOD=THE  RATIO  BETWEEN  THE  ABOVE  d  AND  THE  MAIN  REFLECTOR 

41  C  DIAMETER 

42  C  ALPHA=THE  ANGLE  BETWEEN  THE  LOWER  EDGE  OF  THE  SUBREFLECTOR 

43  C  TO  THE  SECOND  FOCUS  (NOT  THE  ORIGIN) 

44  C  BETA=THE  ANGLE  BETWEEN  THE  AXIS  OF  SYMMETRY  AND  THE 

45  C  IMAGINARY  LINE  BETWEEN  THE  ORIGIN  TO  THE  EDGE  OF  THE 

46  C  SUBREFLECTOR 

47  C  ECC=ECCENTRICITY  OF  THE  SUBREFLECTOR  (ELLIPSOID) 

48  C 

49  C 

50  C  FIRST  GENERATE  PARABOLOID  CONTOUR 

51  C 
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1 

2 

3 

4 

5 

6 

7 

8 
9 

10 

11 

12 

13 

14 

15 

16 

17 

18 

19 

20 
21 
22 

23 

24 

25 

26 

27 

28 

29 

30 

31 

32 

33 

34 

35 

36 

37 

38 

39 

40 

41 

42 

43 

44 

45 

46 

47 

48 

49 

50 

51 


WRITE (6,*)  'ENTER  ALPHA' 

READ (  5  ,  * )  ALPHA 
WRITE (6,*)  'ENTER  BETA' 

READ ( 5 , * )  BETA 

WRITE (6,*)  'ENTER  DM  (IN  WAVELENGTH)' 

READ ( 5 , * )  DM 
WRITE (6,*)  'FOD' 

READ (5,*)  FOD 
WRITE (6,*)  'dOD' 

READ ( 5 , * )  dOD 

ALPHA=ALPHA*RAD 

BETA=BETA*RAD 

FM=FOD*DM 

d=dOD*DM 

Rmax=FM-d/2 ./TAN (PI-ALPHA) 

WRITE (6,*)  'MAXIMUM  FEEDER  DISTANCE  IS',  Rmax 
WRITE (6,*)  'ENTER  Rfeed  (IN  WAVELENGTH)' 

READ (5,*)  Rfeed 
C 

C  FM=FOCUS  OF  THE  MAIN  REFLECTOR  (PARABOLOID) 

C  ZSHIFT=INITIALIZE  THE  ORIGIN  TO  BE  AT  ZERO 
C 

ZSHIFT=FM-Rf eed 
DO  5  1=1, NP1 

TH=FLOAT ( I - 1 ) *ALPHA/FLOAT (NP1-1 ) 

RM=2 . *FM/ ( 1 . +COS (TH) ) 

ZH(I)=(ZSHIFT-RM*COS(TH) )*BK 
RH ( I ) = (d/2 . +RM*SIN (TH) ) *BK 
5  CONTINUE 

C 

C  SUBREFLECTOR  CONTOUR 

C 

PHIZ=ATAN (d/2 ./ZSHIFT) 

DELTA=PI -ALPHA-PHI Z 
CPZ=COS (PHIZ) 

CPS=COS (BETA) 

BB=2*SIN (DELTA) /SIN (ALPHA) 

AA=bb*COS (PHIZ) -1 
ECC= (BB-SQRT (BB**2-4* (AA) ) )/2/AA 
WRITE (6,*)  'ECCENTRICITY  =' , ECC 
EDGEE=CPS** ( FEXP*2 . ) 

EDGEH=CPS** (FHXP*2 . ) 

EEDB=10 . *ALOG10 (EDGEE) 

EHDB=10 . *ALOG10 (EDGEH) 

PSPILL=FH2* ( 1 . -EDGEE) +FE2  * ( 1 . -EDGEH) 
PSPILL=PSPILL/ (FH2+FE2) 

PSDB=10 . *ALOG10 ( PS PILL) 

DO  6  1=1, NP2 

THETA=FLOAT ( I — 1 ) * (BETA) /FLOAT (NP2-1 ) 

RSS=AA* ( (1. -ECC**2. )/ (l.-ECC*COS (PHIZ-THETA) ) ) *BK 
ZH ( I+NP1 ) =RSS*COS (THETA) 
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26 
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28 

29 
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RH ( I+NP1) =RSS*SIN (THETA) 

6  CONTINUE 

OPEN  (8  ,  f  ile=  '/home/srvl/vered/thesis/outADE  ' ) 

WRITE (8,8000)  DM, RS , ECC , PHI Z/RAD , ALPHA/RAD , NP1 , NP2 , FEXP , 
*  FHXP 

8000  FORMAT (//, 5X, ' ***  ADE  ANTENNA  SYSTEM  ***',//, 5X, 

*' REFLECTOR  PARAMETERS  (LENGTHS  IN  WAVELENGTHS 5X , 
* • MAIN  REFL  DIA= ' , F8 . 3 , / , 5X , ' SUB  REFL  DISTANCE= ' , F8 . 3 , / , 
*5X, ' ECCENTRICITY= ' ,F8.4,/,5X, 'PHIZ  ( DEG ) = ' ,F8.2,/,5X, 
*' ALPHA  ( DEG)  =  ' ,F8.2,/,5X,  'MAIN  REFL  POINTS  (NP1)=' ,14,/, 
*5X , ' SUB  REFL  POINTS  (NP2 )=', 14 ,/, 5X , ' E-PLANE  FEED 
*EXPONENT= ' , F8 . 3 , / , 5X , ' H-PLANE  FEED  EXPONENT= ' , F8 . 3 ) 
WRITE (8,30)  NT, NPHI 

30  FORMAT (/,7X, '  NT  NPHI ' , / , 5X , 15 , 2X , 15) 

IF(ISEG.EQ.O)  WRITE(8, 1300) 

1300  FORMAT (/,5X, 'INDEX' ,8X, 'Z(I) ' , 10X, 'RHO(I) ' ) 

OPEN ( 1 , f ile= ' /home/srvl/vered/matlab/zeem.m ' ) 

OPEN (2 , f ile= '/home/srvl/vered/matlab/rhom.m' ) 

OPEN(3, file='/home/srvl/vered/matlab/zees.m' ) 

OPEN ( 4 , f ile= ' /home/srvl/vered/matlab/rhos . m ' ) 

DO  52  1=1, NP 

IF(ABS(ZH(I) ) .LT. .001)  ZH(I)=0. 

IF(ABS (RH(I) ) .LT. . 001)  RH(I)=0. 

ZHB=ZH(I)/BK 
RHB=RH ( I ) / BK 
IF(I.LE.NPl)  THEN 
WRITE (1,31)  ZHB 
WRITE (2,31)  RHB 
ENDIF 

IF ( I . GE . NP1+1 )  THEN 
WRITE ( 3 , 31)  ZHB 
WRITE (4,31)  RHB 
ENDIF 

51  if (iseg.eq.O)  WRITE ( 8 , 8004 )  I, ZHB, RHB 

52  CONTINUE 

8004  FORMAT (6X,I4,4X,F8.3,8X,F8.3) 

31  FORMAT (F8. 4) 

CLOSE ( 1 ) 

CLOSE ( 2 ) 

CLOSE (3) 

CLOSE ( 4 ) 

write (6,*)  'geometry  defined' 

I F ( I CALC . EQ . 1 )  THEN 

open ( 3 , f ile= '/home/srvl/vered/thesis/current ' ) 

READ ( 3 , * )  DM0 , DSO , FCO , PHIVO , PHIRO 
READ ( 3 , * )  ECCO , FODO , NP10 , NP2 0 , FEXPO , FHXPO 
WRITE ( 6 , * ) DM0 , DSO , FCO , PHIVO , PHIRO , ECCO , FODO , NP10 , NP20 , 
*FEXF0 

WRITE (6, *)  DM, DS , FC , PHIV, PHIR , ECC , FOD, NP1 , NP2 , FEXP 
READ ( 3 , * )  (RH ( J ) , J=1 , NP) 

READ ( 3 , * )  ( ZH ( J) , J=1 , NP) 
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READ ( 3 , * )  (C ( J) , J=1 ,  N) 

WRITE (6,*)  'CURRENTS  READ  SUCCESSFULLY' 

CLOSE (3) 

ENDIF 

I F ( I CALC . EQ . 0 )  THEN 

CALL  SPHERE ( FEXP , FHXP , NP1 , NP2 , NT , RH , ZH , XT , AT , BK , R) 

DO  22  J=1,N 
B(J)=R(J) 

IF(J.GT.MT)  B(J)=-B(J) 

22  CONTINUE 

write (6,*)  'excitation  coefficients  computed' 

ENDIF 

CALL  ZMAT (1,1, NP1 , NP2 , NPHI , NT , RH , ZH , X , A , XT , AT , Z ) 
write (6,*)  'back  from  zmat' 

C 

C  NULL  OUT  ROWS  AND  COLUMNS  CORRESPONDING  TO  DUMMY  SEGMENTS 
C 

DO  700  LSS  =  1 , N 
LS=LSS-1 

Z (LS*N+NP1)=(0. ,0.) 

Z (LS*N+NP1-1) = (0 • , 0. ) 

700  Z (LS*N+NP1+MT) =(0. , 0. ) 

DO  701  LS=1 , N 

Z ( (NP1-2) *N+LS)=(0. ,0.) 

Z ( (NP1-1) *N+LS ) = ( 0 . , 0. ) 

701  Z( (NP1-1+MT) *N+LS)=(0. ,0. ) 

Z ( (NP1-2) *N+NP1-1)=(1. ,0. ) 

Z ( (NP1-1) *N+NP1)=(1. , 0. ) 

Z( (MT+NP1-1) *N+NP1+MT)=(1. ,0.) 

CALL  DECOMP (N, IPS, Z) 

CALL  SOLVE (N, IPS, Z,B,C) 

write(6,*)  'coefficients  calculated' 

IF(ICWRT.EQ.O)  THEN 

OPEN ( 3 , f ile= '/home/srvl/vered/thesis/current ' ) 

WRITE ( 3 , * )  DM , DS , FC , PHIV , PHIR , ECC 
WRITE (3,*)  FOD,NPl,NP2, FEXP, FHXP 
WRITE ( 3 , * )  (RH ( J) , J=1 , NP) 

WRITE ( 3 , * )  ( ZH ( J ) , J=1 , NP) 

WRITE ( 3 , * )  (C ( J) , J=1 , N) 

CLOSE (3) 

ENDIF 

991  CONTINUE 

IT=INT ( (STOP-START) /DT) +1 
C 

C  RECEIVER  PHI  CUTS:  DO  0  AND  90  (E-  AND  H-  PLANES) 

C 

DO  400  IP=1 , 2 
ECX=0 . 

PH 1=0 . 

IF ( IP . EQ . 2 )  PHI=90 . 

PHR=PHI*RAD 
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CNP=COS ( PHR) 

SNP=SIN (PHR) 

DO  500  1=1 , IT 
THETA=FLOAT ( I - 1 ) *DT+START 
THX=THETA*RAD 

CALL  PLANE ( 1 , 1 , NP1 , NP2 , NT , RH , ZH , XT , AT , THX , R) 

ETST= ( 0 . ,0. ) 

EPST= ( 0 . ,0. ) 

ETSP= ( 0 . ,0. ) 

EPSP= ( 0 .  ,0.  ) 

DO  100  10=1, MT 
ETST=ETST+R(I0) *C(I0) 

100  ETSP=ETSP+R ( IO+N) *C (10) 

DO  200  10=1, MP 
EPST=EPST+R ( IO+MT) *C(I0+MT) 

200  EPSP=EPSP+R ( IO+N+MT) *C(I0+MT) 

ETF=CNP*FE ( FEXP , THX ) 

EPF=SNP*FH ( FHXP , THX ) 

ETT=2 . * (ETST+EPST) *UC*CNP 

EPT=2 . *U* (ETSP+EPSP) *UC*SNP 

ET=ETF+ETT 

EP=EPF+EPT 

EC= (ET*CNP-EP*SNP) 

EX= (ET*SNP+EP*CNP) 

ECV ( I ) =CABS (EC) 

EXV ( I ) =CABS ( EX ) 

ECR=REAL ( EC ) 

ECI=AIMAG (EC) 

EXR=REAL ( EX ) 

EXI=AIMAG (EX) 

PHC ( I ) =ATAN2 (ECI , ECR+1 . e-10) /RAD 
PHX ( I ) =ATAN2 (EXI , EXR+1 .e-10) /RAD 
ANG ( I ) =THETA 

ECX=AMAX1 (ECX , ECV ( I ) ,EXV(I) ) 

500  CONTINUE 

IF(IP.EQ.l)  THEN 

WRITE (8, 103)  ECX, EEDB, EHDB, PSPILL, PSDB 
103  FORMAT (/,5X, 'E-PLANE  EDGE  TAPER= ' , F8 . 2 , 

*  /,5X, 'H-PLANE  EDGE  TAPER= ' , F8 . 2 , / , 5X , ' TOTAL  RAD 

*  POWER= ' , E15 . 5 , 

*  /, 5X, 'SPILLOVER  LOSS  = ' , F8 . 4 , / , 5X, '  IN  dB 

*  =' , f8.2) 

ENDIF 

WRITE (8,107)  PHI 

107  FORMAT (/, 5X, 'PHI  OF  RECEIVER  ( DEG) = ' , F8 . 2 ) 

DO  600  1=1, IT 

ECP ( I )  = (ECV ( I ) /ECX) **2 

EXP ( I ) = (EXV ( I ) /ECX) **2 

ECP ( I ) =AMAX1 ( ECP ( I ) , .00001) 

EXP ( I ) =AMAX1 (EXP ( I ) , .00001) 

ECP ( I ) =10 . *ALOG10 (ECP ( I ) ) 


52 


1  EXP (I) =10 . *ALOG10 (EXP (I) ) 

2  600  CONTINUE 

3  IF ( I PRINT . NE . 0)  GO  TO  300 

4  WRITE(8, 5015) 

5  5015  FORMAT (/,7X, 'ANGLE' , 15X, 'CO-POL' ,25X, 'X-POL' ,/,7X, 

6  1' (DEG) * ,4X, 2 ( ' (VOLTS) ' ,  4X, ' (DEG) ' ,3X, * (DB-REL) ' ,  4X) ) 

7  DO  9000  L=1 , IT 

8  WRITE(8, 5016) ANG(L) ,ECV(L) ,PHC(L) ,ECP(L) ,EXV(L) ,PHX(L) 

9  1 , EXP ( L) 

10  5016  FORMAT (5X,F6.1,3X,2(F8.2,3X,F7.1,3X,F7.2,3X)) 

11  9000  CONTINUE 

12  300  CONTINUE 

13  IF (IP. EQ. 1)  THEN 

14  OPEN  (2  ,  file='/home/srvl/vered/matlab/cang.m' ) 

15  OPEN(3, file='/home/srvl/vered/matlab/ccpole.m' ) 

16  OPEN (4 , f ile= '/home/srvl/vered/matlab/cxpole . m ' ) 

17  DO  9097  1=1, IT 

18  WRITE (2,5019)  ANG(I) 

19  WRITE (3 , 5019)  ECP(I) 

20  9097  WRITE(4 , 5019)  EXP(I) 

21  CLOSE (2) 

22  CLOSE (3) 

23  CLOSE (4) 

24  ENDIF 

25  IF ( IP. EQ . 2 )  THEN 

26  OPEN ( 3  ,  file^/home/srvl/vered/matlab/ccpolh.m' ) 

27  OPEN ( 4 , file= '/home/srvl/vered/matlab/cxpolh.m' ) 

28  DO  9098  1=1, IT 

29  WRITE (3 , 5019)  ECP(I) 

30  9098  WRITE(4 , 5019)  EXP(I) 

31  CLOSE (3) 

32  CLOSE ( 4 ) 

33  ENDIF 

34  400  CONTINUE 

35  5019  FORMAT ( F8 . 3 ) 

36  C 

37  C  program  gain.f 

38  C 

39  C  FAR  FIELD  PATTERN  INTEGRATION  FOR  LOSS  CALCULATIONS. 

40  C  USES  THE  PREVIOUSLY  CALCULATED  CURRENTS  TRANSFERRED  FROM 

41  C  pardip. f . 

42  C  ***  general  wire/bor  geometries  can  be  handled  *** 

43  C  >>>>  INPUT  FILE  IS  current 

44  C  >»>  HAS  A  GAIN  NORMALIZATION  LOOP  INCORPORATED 

45  C 

46  OPEN ( 1 , FILE= ' /home/srvl/vered/thesis/outgaus ' ) 

47  READ ( 1 , * )  NPHI 

48  NT2=NT/2 

49  DO  7  K=1 , NT2 

50  K1=NT-K+ 1 

51  AT(K1) =AT(K) 
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XT (Kl) =XT (K) 

7  XT ( K) =-XT ( K) 

DO  8  K=1 , NPHI 

8  READ ( 1 ,  * )  XG(K),AG(K) 

NDIVT=1 

NDIVP=1 

S1=0. 

S2=180 . 

Q1=0 . 

Q2=180. 

C 

C  GAIN  NORMALIZATION  LOOP  (USE  THETA=PHI=0  AND  THETA  POL) 

C 

CNP=1 . 

CALL  PLAN (1,1, NP , NT , RH , ZH , XT , AT , 0 . ,R) 

ETST= ( 0 . , 0 . ) 

EPST= ( 0 . ,0. ) 

DO  101  10=1, MT 
101  ETST=ETST+R ( 10) *C(I0) 

DO  201  10=1 , MP 

201  EPST=EPST+R ( IO+MT) *C ( IO+MT) 

ETT=2 . * (ETST+EPST) *UC*CNP 
ETF=CNP*FE (FEXP, 0. ) 

ET=ETT+ETF 

EMAX=CABS(ET) 

WRITE (6,*)  ' EMAX= ' , EMAX 

WRITE (8,8008 )  NT, NPHI , NDIVP, NDIVT , SI , S2 , Q1 , Q2 
8008  FORMAT (//, 10X, ' NT= ' ,I3,/,10X, 'NPHI=' ,I4,/,10X, 

*  'NUMBER  OF  INTEGRATION  INTERVALS  IN  PHI= ' , i 3 , / , lOx , 

*  'NUMBER  OF  INTEGRATION  INTERVALS  IN  THETA= ' , i3 , / , lOx , 

*  'THETA  LIMITS:  START , STOP= ' , F7 . 2 , 2X , F7 . 2 , / , 10X , 

*  'PHI  LIMITS:  START , STOP= ' , F7 . 2 , 2X , F7 . 2 ) 
c 

c  theta  integration  intervals 
c 

S1=S1*RAD 
S2=S2  *RAD 
Q1=Q1*RAD 
Q2=Q2*RAD 

DS=(S2-S1) /FLOAT (NDIVT) 

DO  502  1=1 , NDIVT+ 1 
S(I)  =  (I-1)  * DS 

WRITE (6,*)  "I,S(I)=", I,S (I) 

502  CONTINUE 
C 

c  phi  integration  intervals 
c 

DQ= (Q2-Q1) /FLOAT (NDIVP) 

DO  503  1=1 , NDIVP+ 1 
Q(I)  =  (I-1)*DQ 

WRITE ( 6 , * )  "I,Q(I)=",I,Q(I) 
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1 

503 

CONTINUE 

2 

SUM=0. 

3 

do  2000  JJ=1 , NDI VP 

4 

Pl=DQ/2 . 

5 

P2=(Q(JJ+1) +Q(JJ) )/2. 

6 

DO  2000  J=1 , NPHI 

7 

PHR=pl*xg(j)+p2 

8 

DO  1000  11=1, NDIVT 

9 

Tl=DS/2 . 

10 

T2=(S(II+1)+S(II) )/2. 

11 

DO  1000  1=1, NPHI 

12 

THR=T1*XG ( I ) +T2 

13 

IF (THR. LT. PI)  GO  TO  510 

14 

THR= ( BK-THR ) 

15 

PHR=PHR+PI 

16 

510 

CONTINUE 

17 

SNP=SIN (PHR) 

18 

CNP=COS (PHR) 

19 

CALL  PLAN (1,1, NP, NT , RH , ZH , XT , AT , THR, R) 

20 

STRA=SIN (THR) *AG(I) *AG(J) 

21 

ETST= ( 0 . , 0. ) 

22 

EPST= ( 0 . , 0. ) 

23 

ETSP= ( 0 . , 0. ) 

24 

EPSP= ( 0 .  ,  0.  ) 

25 

DO  104  10=1 ,MT 

26 

ETST=ETST+R(IO) *C(I0) 

27 

104 

ETSP=ETSP+R ( IO+N) *C ( 10) 

28 

DO  204  10=1, MP 

29 

EPST=EPST+R ( IO+MT) *C( IO+MT) 

30 

204 

EPSP=EPSP+R ( IO+N+MT) *C(I0+MT) 

31 

ETT=2 . * (ETST+EPST) *UC*CNP 

32 

EPT=2 . *U* (ETSP+EPSP) *UC*SNP 

33 

ETF=CNP*FE (FEXP, THR) 

3  4 

EPF=SNP*FH ( FHXP, THR) 

35 

ET=ETT+ETF 

36 

EP=EPT+EPF 

37 

SUM=SUM+ (CABS (ET) **2+CABS (EP) **2) *STRA 

38 

1000 

CONTINUE 

39 

2000 

CONTINUE 

40 

PRAD=T1*P1*SUM 

41 

GAIN=2 . *PI*EMAX**2/PRAD 

42 

GDB=10 . *aloglO (GAIN) 

43 

EFF=GAIN/ (PI*DM) **2 .  *100. 

44 

WRITE (8, 8010)  PRAD , EMAX , GAIN , GDB , EFF 

45 

8010 

FORMAT (//, 1 OX, ’TOTAL  RADIATED  POWER= ’ , E15 . 8 ,/ , 10X, 

46 

*  'MAXIMUM  FIELD  VALUE  ( V/M) = ' , E15 . 8 , / , 10X , 

47 

*  ' GAIN= ' , E15 . 8 , / , lOx , ' (dB) = ' , f 8 . 2 , 

48 

*  1 EFFICIENCY= ',E15.8,2x,  '%') 

49 

WRITE ( 6 , * ) "PRAD, GAIN , dB= , EFFICIENCY" , PRAD, GAIN, GDB, EFF 

50 

32 

FORMAT ( F8 . 4 ) 

51 

9999 

STOP 
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END 

SUBROUTINE  ZMAT (Ml , M2 , NP1 , NP2 , NPHI , NT , RH , ZH , X , A , XT , AT , Z ) 
C 

C  subroutine  to  compute  the  impedance  elements  for  BORs. 

C  Original  version  from  Mautz  and  Harrington  (Improved  E-Field 
C  Method  . . . ) 

C  THE  SUBROUTINE  ZMAT  CALLS  THE  FUNCTION  BLOG 
C 

COMPLEX  U1/U2>U3>U4,U5>U6,U7,U8, U9 , UA, UB, G4A (10) ,G5A(10) 
COMPLEXG  6 A ( 1 0 ) ,G4B(10) ,G5B(10) ,G6B(10) , H4A, H5A, H6A , H4B 
COMI  ..EX  CMPLX ,  H6B ,  UC ,  UD ,  GA  (100)  ,GB(100)  ,H5B,  Z  (300000) 
DIMENSION  RH (300) ,ZH(300) ,X(100) ,A(100) ,XT(10) ,AT(10) 
DIMENSION  RS (300) ,ZS(300) ,D(300) ,DR(300) ,DZ(300) 
DIMENSION  DM  (3  00)  ..  C2  (100)  , C3 (100)  ,  R2  ( 10 )  ,  Z2  ( 10 )  ,  R8  ( 10 ) 
DIMENSION  C4 ( 100 ) ,C5(100) ,C6(100) f Z7 ( 10)  , R7 (10)  , Z8  ( 10) 
NP=NP1+NP2 
CT=2  . 

CP=.  1 

DO  10  1=2 , NP 
12=1-1 

RS ( 12 ) = . 5* (RH ( I ) +RH ( 12 ) ) 

ZS(I2)=.5*(ZH(I)+ZH(I2) ) 

Dl= . 5* (RH ( I ) -RH ( 12 ) ) 

D2= . 5* ( ZH ( I ) -ZH ( 12 ) ) 

D ( 12 ) =SQRT (D1*D1+D2*D2) 

DR(I2)=D1 
DZ ( 12 ) =D2 

I F ( RS (12). EQ . 0 . )  RS ( 12 ) =1 . 

DM ( 12 ) =D ( 12 ) /RS ( 12 ) 

10  CONTINUE 
M3=M2-M1+1 
M4=M1-1 
PI2=1. 570796 
DO  11  K=1 , NPHI 
PH=PI2* (X(K)+1.) 

C2 (K) =PH*PH 
SN=SIN( . 5* PH) 

C3 (K) =4 . *SN*SN 
A1=PI2*A(K) 

D4=.5*A1*C3 (K) 

D5=Al*COS (PH) 

D6=A1*SIN ( PH) 

M5=K 

DO  29  M=1 , M3 
PHM= (M4+M) *PH 
A2=COS ( PHM) 

C4 (M5)=D4*A2 
C5(M5)=D5*A2 
C6 (M5 )=D6*SIN( PHM) 

M5=M5+NPHI 
29  CONTINUE 
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1 

11 

CONTINUE 

2 

MP=NP-1 

3 

MT=MP- 1 

4 

N=MT+MP 

5 

N2N=MT*N 

6 

N2=N*N 

7 

U1=(0. , .5) 

8 

U2= ( 0 . ,2.) 

9 

JN=-1-N 

10 

DO  15  JQ=1,MP 

11 

KQ=2 

12 

IF(JQ.EQ.l)  KQ=1 

13 

IF(JQ.EQ.MP)  KQ=3 

14 

R1=RS (JQ) 

15 

Z1=ZS ( JQ) 

16 

D1=D ( JQ) 

17 

D2=DR ( JQ) 

18 

D3=DZ (JQ) 

19 

D4=D2/R1 

20 

D5=DM ( JQ) 

21 

SV=D2/D1 

22 

CV=D3/D1 

23 

T6=CT*D1 

24 

T62=T6+D1 

25 

T62=T62  *T62 

26 

R6=CP*R1 

27 

R62=R6*R6 

28 

DO  12  L=1 , NT 

29 

R2 ( L) =R1  +  D2  *XT ( L) 

30 

Z2 ( L) =Z 1  +  D3  *XT ( L) 

31 

12 

CONTINUE 

32 

U3=D2  *U1 

33 

U4  =  D3  *U1 

34 

DO  16  IP=1,MP 

35 

R3=RS (IP) 

36 

Z3=ZS (IP) 

37 

R4=R1-R3 

38 

Z4=Z1-Z3 

39 

FM=R4  *SV+Z4  *CV 

40 

PHM=ABS (FM) 

41 

PH=ABS (R4*CV-Z4*SV) 

42 

D6=PH 

43 

IF(PKM.LE.Dl)  GO  TO  26 

44 

D6=PHM-D1 

45 

D6=SQRT (D6*D6+PH*PH) 

46 

26 

IF ( IP. EQ . JQ)  GO  TO  27 

47 

KP=1 

48 

I F (T6 . GT . D6 )  KP=2 

49 

IF (R6 . GT. D6)  KP=3 

50 

GO  TO  28 

51 

27 

KP=4 
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1 

28 

GO  TO  (41,42,41,42) ,KP 

2 

42 

DO  40  L=1 , NT 

3 

D7=R2 (L) -R3 

4 

D8=Z2 (L) -Z3 

5 

Z7 (L)  =D7*D7  +  D8  *D8 

6 

R7 ( L) =R3  *R2 (L) 

7 

Z8 (L) =. 25*Z7 (L) 

8 

R8(L)=.25*R7(L) 

9 

40 

CONTINUE 

10 

Z4=R4*R4+Z4*Z4 

11 

R4=R3  *R1 

12 

R5= . 5*R3  *SV 

13 

DO  33  K=1 , NPHI 

14 

A1=C3 (K) 

15 

RR=Z4+R4  *A1 

16 

UA=0 . 

17 

UB=0. 

18 

IF (RR . LT . T62 )  GO  TO  34 

19 

DO  35  L=1,NT 

20 

R=SQRT(Z7 ( L) +R7 ( L) *A1) 

21 

SN=-SIN (R) 

22 

CS=COS (R) 

23 

UC=AT ( L) /R*CMPLX ( CS , SN ) 

24 

UA=UA+UC 

25 

UB=XT(L) *UC+UB 

26 

35 

CONTINUE 

27 

GO  TO  36 

28 

34 

DO  37  L=1 , NT 

29 

R=SQRT(Z8 (L) +R8 (L) *A1) 

30 

SN=-SIN(R) 

31 

CS=COS (R) 

32 

UC=AT ( L) /R*SN*CMPLX ( -SN , CS ) 

33 

UA=UA+UC 

34 

UB=XT ( L) *UC+UB 

35 

37 

CONTINUE 

36 

A2=FM+R5*A1 

37 

D9=RR-A2  *A2 

38 

R=ABS ( A2 ) 

39 

D7=R-D1 

40 

D8=R+D1 

41 

D6=SQRT (D8*D8+D9 ) 

42 

R=SQRT(D7*D7+D9) 

43 

IF ( D7 . GE . 0 . )  GO  TO  38 

44 

Al=ALOG(  (D8+D6)  *  ( -D7-f-R)  /  D9  ) /D1 

45 

GO  TO  39 

46 

38 

Al=ALOG( (D8+D6 ) / (D7+R) )/Dl 

47 

39 

UA=A1+UA 

48 

UB=A2  * ( 4 . / ( D6+R) -A1)/D1+UB 

49 

36 

GA (K) =UA 

50 

GB (K) =UB 

51 

33 

CONTINUE 
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1 

K1=0 

2 

DO  45  M=1,M3 

3 

H4A=0 . 

4 

H5A=0 . 

5 

H6A=0. 

6 

H4B=0 . 

7 

H5B=0 . 

8 

H6B=0. 

9 

DO  46  K=1 , NPHI 

10 

K1=K1+1 

11 

D6=C4 (Kl) 

12 

D7=C5 (Kl) 

13 

D8=C6 (Kl) 

14 

UA=GA (K) 

15 

UB=GB(K) 

16 

H4A=D6*UA+H4A 

17 

H5A=D7*UA+H5A 

18 

H6A=D8*UA+H6A 

19 

H4B=D6*UB+H4B 

20 

H5B=D7*UB+H5B 

21 

H6B=D8*UB+H6B 

22 

46 

CONTINUE 

23 

G4A (M) =H4 A 

24 

G5A (M) =H5A 

25 

G6A (M) =H6A 

26 

G4B (M) =H4B 

27 

G5B (M) =H5B 

28 

G6B (M) =H6B 

29 

45 

CONTINUE 

30 

IF (KP. NE . 4 )  GO  TO  47 

31 

A2=D1/ (PI2*R1) 

32 

D6=2 ./Dl 

33 

D8  =  0 . 

34 

DO  63  K=1 , NPHI 

35 

A1=R4  *C2 (K) 

36 

R=R4  *C3 (K) 

37 

IF (R. LT. T62 )  GO  TO  64 

38 

D7  =  0 . 

39 

DO  65  L=1 , NT 

40 

D7=D7+AT ( L) /SQRT (Z7 ( L)  +A1) 

41 

65 

CONTINUE 

42 

GO  TO  66 

43 

64 

A1=A2/ (X (K) +1 . ) 

44 

D7=D6*ALOG (Al+SQRT ( 1 . +A1*A1) ) 

45 

66 

D8=D8+A(K) *D7 

46 

63 

CONTINUE 

47 

Al= . 5*A2 

48 

A2-1./A1 

49 

D'i--PI2*D8  +  2  ./Rl*  (BLOG(A2)  +A2 

50 

DO  67  M=1 , M3 

51 

G5A (M) =D8+G5A (M) 
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1 

67 

CONTINUE 

2 

GO  TO  47 

3 

41 

DO  25  M=1 , M3 

4 

G4 A (M) =0 . 

5 

G5A (M) =0  . 

6 

G6A (M) =0. 

7 

G4B (M) =0 . 

8 

G5B (M) =0. 

9 

G6B (M) =0. 

10 

25 

CONTINUE 

11 

DO  13  L=1 , NT 

12 

A1=R2 (L) 

13 

R4=A1-R3 

14 

Z4=Z2 (L) -Z3 

15 

Z4=R4*R4+Z4*Z4 

16 

R4=R3  *A1 

17 

DO  17  K=1 , NPHI 

18 

R=SQRT (Z4+R4*C3 (K) ) 

19 

SN=-SIN(R) 

20 

CS=COS (R) 

21 

GA (K) =CMPLX ( CS , SN) /R 

22 

17 

CONTINUE 

23 

D6=0. 

24 

IF (R62 . LE . Z4 )  GO  TO  51 

25 

DO  62  K=1 , NPHI 

26 

D6=D6+A (K) /SQRT ( Z4+R4  *CZ (K) ) 

27 

62 

CONTINUE 

28 

Z4=3 . 14 159 3/ SQRT (Z4/R4 ) 

2  3 

D6=-PI2  *D6+ALOG ( Z4  +  SQRT (l.+Z4*Z4) )/SQRT(R4) 

30 

51 

A1=AT  ( L) 

31 

A2=XT(L) *A1 

32 

K1=0 

33 

DO  30  M=1,M3 

34 

U5=0 . 

35 

U6=0 . 

36 

U7=0 . 

37 

DO  32  K=1 , NPHI 

38 

UA=GA(K) 

39 

K1=K1+1 

40 

U5=C4 (Kl) *UA+U5 

41 

U6=C5(K1) *UA+U6 

42 

U7=C6 (Kl ) *UA+U7 

43 

32 

CONTINUE 

44 

U6=D6+U6 

45 

G4A (M) =A1*U5+G4A (M) 

46 

G5A(M)=A1*U6+G5A(M) 

47 

G6A (M) =A1*U7+G6A (M) 

48 

G4B(M)=A2*U5+G4B(M) 

49 

G5B(M)=A2*U6+G5B(M) 

50 

G6B (M) =A2  *U7+G6B (M) 

51 

30 

CONTINUE 
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29 
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31 
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13  CONTINUE 

47  A1=DR ( IP) 

UA=A1*U3 
UB=DZ (IP) *U4 
A2=D (IP) 

D6  =  -A2  *D2 
D7=D1 *A1 
D8=D1 *A2 
JM=JN 

DO  31  M=1 , M3 
FM=M4+M 
A1=FM*DM (IP) 

H5A=G5A (M) 

H5B=G5B(M) 

H4A=G4A (M) +H5A 
H4B=G4B(M) +H5B 
H6A=G6A (M) 

H6B=G6B (M) 

U7=UA*H5A+UB*H4 A 

U8=UA*H5B+UB*H4  B 

U5=U7-U8 

U6=U7+U8 

U7=-U1*H4A 

U8=D6*H6A 

U9=D6*H6B-A1*H4A 

UC=D7* (H6A+D4  *H6B) 

UD=FM*D5*H4A 

K1=IP+JM 

K2=K1+ 1 

K3=K1+N 

K4=K2+N 

K5=K2+MT 

K6=K4+MT 

K7=K3+N2N 

K8=K4+N2N 

GO  TO  (18,20,19) , KQ 

18  Z(K6)=U8+U9 

IF (IP. EQ. 1)  GO  TO  21 
Z ( K3 ) =Z ( K3 ) +U6-U7 
Z (K7) =Z (K7 ) +UC-UD 
IF (IP. EQ.MP)  GO  TO  22 

21  Z (K4 ) -U6+U7 
Z (K8)=UC+UD 
GO  TO  22 

19  Z ( K5 ) =Z (K5 ) +U8-U9 

IF ( IP. EQ . 1)  GO  TO  23 
Z (Kl) =Z (Kl) +U5+U7 
Z (K7 ) =Z ( K7 ) +UC-UD 
IF(IP. EQ.MP)  GO  TO  22 

23  Z ( K2 ) =Z (K2 ) +U5-U7 
Z (K8)=UC+UD 
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1 

GO  TO  22 

2 

20 

Z  (K5)=Z (K5) +U8-U9 

3 

Z (K6) =U8+U9 

4 

IF(IP.EQ.l)  GO  TO  24 

5 

Z (K1)=Z (Kl) +U5+U7 

6 

Z  (K3 )  =Z  (K3 )  -+U6-U7 

7 

Z (K7) =Z (K7) +UC-UD 

8 

IF(IP.EQ.MP)  GO  TO  22 

9 

24 

Z (K2 ) =Z (K2 ) +U5-U7 

10 

Z (K4)=U6+U7 

11 

Z(K8)=UC+UD 

12 

22 

Z (K8+MT) =U2  * (D8  * (H5A+D4*H5B) -A1*UD) 

13 

JM=JM+N2 

14 

31 

CONTINUE 

15 

16 

CONTINUE 

16 

JN=JN+N 

17 

15 

CONTINUE 

18 

DO  100  LSS=1 , N 

19 

LS=LSS-1 

20 

Z ( LS*N+NP1 )  =  ( 0 .  ,  0.  ) 

21 

Z ( LS*N+NP1-1 )  =  ( 0 .  ,0.  ) 

22 

100 

Z(LS*N+NP1+MT)=(0. ,0. ) 

23 

DO  101  LS  =  1 , N 

24 

Z ( (NP1-2) *N+LS)=(0. ,0.) 

25 

Z ( (NP1-1) *N+LS )  =  ( 0 . , 0.  ) 

26 

101 

Z ( (NP1-1+MT) *N+LS )  =  ( 0 . ,0.  ) 

27 

Z( (NP1-2) *N+NP1-1)=(1. ,0.  ) 

28 

Z( (NP1-1) *N+NP1)  =  (1.  ,0.  ) 

29 

Z ( (MT+NP1-1) *N+NP1+MT)=(1. , 0.  ) 

30 

RETURN 

31 

END 

32 

SUBROUTINE  SPHERE ( FEXP , FHXP , NP1 , NP2 , NT , RH , ZH , XT , AT , BK , R) 

33 

C 

34 

C 

SPHERICAL  WAVE  EXCITATION  VECTOR  -  SOURCE  ON  Z  AXIS 

35 

C 

(r.=+l,-l  ONLY).  DISTANCE  BETWEEN  ORIGINS  (0  TO  O')  IS  ZF 

36 

C 

37 

COMPLEX  R ( 2000 ) , U , U2 , U3 , CEXP , PSI ( 3 00 ) , FIE , F2E , F3E , F4E , 

38 

COMPLEX  CMPLX , UO , UES , UEC , UH , F1H , F2H 

39 

DIMENSION  TH(10) ,RH(300) ,ZH(300) ,XT(10) ,AT(10) 

40 

NP=NP1+NP2 

41 

MP=NP-1 

42 

MT=MP-1 

43 

N=MT+MP 

44 

U= ( 0 . , 1 . ) 

45 

PI=3 . 141593 

46 

PI2=PI/2 . 

47 

DO  12  IP=1 , MP 

48 

K2  =  IP 

49 

I=IP+1 

5n 

DR= . 5* (RH ( I ) -RH (IP) ) 

51 

DZ=. 5* (ZH(I) -ZH(IP) ) 
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D1=SQRT (DR*DR+DZ*DZ) 

R1=.5*(RH(I)+RH(IP) ) 

IF(ABS(R1) .LT.l.E-5)  Rl=l . 
Z1=.5*(ZH(I)+ZH(IP) ) 

SVP=DR/D1 
CVP=DZ/D1 
DO  13  L=1 , NT 
QS=DR*XT ( L) 

QC=DZ*XT ( L) 

TH(L)=ATAN2 ( (Rl+QS) , (Zl+QC) ) 

RP=SQRT( (Rl+QS) **2+(Zl+QC) **2) 

PSI (L)=CEXP(CMPLX(0. ,-RP) ) /RP 
13  CONTINUE 
F1E= ( 0 . ,0. ) 

F2E= ( 0 . ,0. ) 

F3E= ( 0 . ,0. ) 

F4E= ( 0 . ,0. ) 

F1H= ( 0 . ,0. ) 

F2H= ( 0 . , 0. ) 

DO  15  L=1 , NT 
U0=AT ( L) *PSI (L) 

UES=FE ( FEXP , TH ( L) ) *SIN(TH(L)  )  *U0 
UEC=FE ( FEXP, TH ( L) ) *COS (TH ( L) ) *UC 
UH=U0*FH ( FHXP , TH ( L) ) 

F1E=F1E+UES 
F2E=F2E+UEC 
F3E=F3E+UES*XT ( L) 

F4E=F4E+UEC*XT(L) 

F1H=F1H+UH 
F2H=F2H+UH*XT ( L) 

15  CONTINUE 

U2=-D1*PI2* ( CVP*F1E-SVP*F2E) 

U3=-D1*PI2* (CVP*F3E-SVP*F4E) 

K1=K2-1 

R (K2+MT) =D1*U*PI* (F1H+D1*SVP*F2H/R1) 

IF ( IP. EQ . 1)  GO  TO  21 
R(K1)=R(K1)+U2-U3 
IF(IP.EQ.MP)  GO  TO  27 
21  R(K2)=U2+U3 
27  CONTINUE 
12  CONTINUE 
C 

C  SIMPLIFIED  METHOD  FOR  HANDLING  MULTIPLE  BODIES 
C 

R(NPl-l) = (0 . ,0. ) 

R (NP1 )  =  ( 0 .  ,0.  ) 

R (NP1+MT) = ( 0 . ,0. ) 

RETURN 

END 

SUBROUTINE  SOLVE (N, IPS ,UL, B, X) 

COMPLEX  UL( 300000) ,B(500) ,X(500) , SUM 
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DIMENSION  IPS (500) 

NP1=N+1 

IP=IPS(1) 

X  ( 1)  =B  (IP) 

DO  2  1=2, N 
IP=IPS (I) 

IPB=IP 
IM1=I-1 
SUM= ( 0 . , 0 . ) 

DO  1  J=1 , IM1 
SUM=SUM+UL (IP) *X(J) 

1  IP=IP+N 

2  X(I)=B(IPB) -SUM 
K2=N* (N-l) 

IP=IPS (N) +K2 

X (N) =X (N) /UL (IP) 

DO  4  IBACK=2 , N 
I=NP1-IBACK 
K2=K2 -N 
IPI=IPS ( I ) +K2 
IP1=I+1 
SUM= ( 0 . , 0 . ) 

IP=IPI 

DO  3  J=IP1 , N 
IP=IP+N 

3  SUM=SUM+UL (IP) *X ( J ) 

4  X(I)=(X(I)-SUM)/UL(IPI) 

RETURN 

END 

SUBROUTINE  DECOMP (N , IPS , UL) 

COMPLEX  UL(300000) , PIVOT, EM 
DIMENSION  SCL(500) ,IPS(500) 

DO  5  1=1, N 
IPS (I) =1 
RN=0 . 

J1=I 

DO  2  J=1 , N 

ULM=ABS(REAL(UL(J1) ) ) +ABS ( AIMAG (UL ( J1 ) ) ) 
J1=J1+N 

IF (RN-ULM)  1,2,2 

1  RN=ULM 

2  CONTINUE 
SCL ( I ) =1 . /RN 

5  CONTINUE 
NM1=N- 1 
K2=0 

DO  17  K=  1 ,  NM 1 
BIG=0 . 

DO  11  I=K, N 
IP=IPS (I) 

IPK=IP*-K2 


64 


1 

2 

3 

4 

5 

6 

7 

8 
9 

10 

11 

12 

13 

14 

15 

16 

17 

18 

19 

20 
21 
22 

23 

24 

25 

26 

27 

28 

29 

30 

31 

32 

33 

34 

35 

36 

37 

38 

39 

40 

41 

42 

43 

44 

45 

46 

47 

48 

49 

50 

51 


SIZE= ( ABS ( REAL (UL ( IPK) ) ) +ABS (AIMAG (UL(IPK) ) ) ) *SCL(IP) 
IF (SIZE-BIG)  11,11,10 

10  BIG=SIZE 
IPV=I 

11  CONTINUE 
IF(IPV-K)  14,15,14 

14  J=IPS (K) 

IPS (K) =IPS ( IPV) 

IPS(IPV)=J 

15  KPP=IPS (K) +K2 
PIVOT=UL ( KPP) 

KP1=K+1 

DO  16  I=KP1,N 
KP=KPP 
IP=IPS ( I ) +K2 
EM=-UL (IP) /PIVOT 
18  UL (IP) =-EM 
DO  16  J=KP1 , N 
IP=IP+N 
KP=KP+N 

UL  ( IP)  =UL  (IP)  +EM*UL  ( KP) 

16  CONTINUE 
K2=K2+N 

17  CONTINUE 
RETURN 
END 

FUNCTION  BLOG (X) 

IF (X . GT . .1)  GO  TO  1 
X2=X*X 

BLOG= ( ( . 075*X2-. 1666667) *X2+1. ) *X 
RETURN 

1  BLOG=ALOG (X+SQRT ( 1 . +X*X) ) 

RETURN 

END 

FUNCTION  FE (FEXP, X) 

PI2=3. 14159/2. 

FE=0. 

IF (ABS (X) .LE.PI2)  FE=ABS (COS (X) ) **FEXP 

RETURN 

END 

FUNCTION  FH ( FHXP , X) 

PI2=3. 14159/2. 

FH=0 . 

IF (ABS (X) .LE.PI2)  FH=-ABS (COS (X) ) **FHXP 

RETURN 

END 

SUBROUTINE  PLANE (Ml , M2 , NP1 , NP2 , NT , RH , ZH , XT , AT , THR , R) 

C 

C  Unchanged  version  from  Mautz  and  Harrington  report. 

C 

COMPLEX  R (2000) ,U ,U1 , UA,UB, FA ( 10) ,FB(10) , F2A, F2B, F1A, FIB 
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1  COMPLEX  U2,U3,U4,U5, cmplx 

2  DIMENSION  RH ( 300) , ZH ( 300) , XT ( 10) , AT ( 10) ,R2 ( 10) , Z2 ( 10) 

3  DIMENSION  BJ(5000) 

4  NP=NP1+NP2 

5  MP=NP-1 

6  MT=MP-1 

7  N=MT+MP 

8  N2=2*N 

9  CC=COS (THR) 


10 

SS=SIN (THR) 

11 

11 

CONTINUE 

12 

U= ( 0 . , 1 . ) 

13 

Ul=3 . 141593*U**M1 

14 

M3=M1+1 

15 

M4=M2+3 

16 

IF(Ml.EQ.O)  M3=2 

17 

M5=Ml+2 

18 

M6=M2+2 

19 

DO  12  IP=1 , MP 

20 

K2  =  IP 

21 

I=IP+1 

22 

DR= . 5* (RH ( I ) -RH (IP) ) 

23 

DZ=. 5* (ZH(I) -ZH(IP) ) 

24 

D1=SQRT (DR*DR+DZ*DZ ) 

25 

Rl= . 25* (RH ( I ) +RH (IP) ) 

26 

IF (R1 . EQ . 0  - )  Rl  =  l . 

27 

Z1=.5*(ZH(I)+ZH(IP) ) 

28 

DR= . 5  *  DR 

29 

D2=DR/R1 

30 

DO  13  L=1 , NT 

31 

R2 (L)=R1+DR*XT(L) 

32 

Z2 (L)=Z1+DZ*XT(L) 

33 

13 

CONTINUE 

34 

D3=DR*CC 

35 

D4=-DZ  *SS 

36 

D5=D1*CC 

37 

DO  23  M=M3 , M4 

38 

FA (M) =0 . 

39 

FB (M) =0 . 

40 

23 

CONTINUE 

41 

DO  15  L=1 , NT 

42 

X=SS*R2 ( L) 

43 

IF (X. GT. . 5E-7 )  GO  TO  19 

44 

DO  20  M=M3 , M4 

45 

BJ (M) =0 . 

46 

20 

CONTINUE 

47 

BJ ( 2 ) =1 . 

48 

S=l. 

49 

GO  TO  18 

50 

19 

M=2 . 8*X+14 . -2 . /X 

51 

I F ( X . LT . . 5 )  M=ll . 8+ALOG10 
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IF (M . LT . M4 )  M=M4 
BJ  (M) =0 . 

JM=M-1 
BJ ( JM) =1 . 

DO  16  J=4 , M 
J2=JM 
JM=JM-1 
J1=JM-1 

BJ ( JM) =J1/X*BJ ( J2 ) -BJ ( JM+2 ) 

16  CONTINUE 
S=0. 

I F ( M . LE . 4 )  GO  TO  24 
DO  17  J=4 , M , 2 
S=S+BJ (J) 

17  CONTINUE 

24  S=BJ ( 2 ) +2 . *S 

18  ARG=Z2 (L)  *CC 

UA=AT ( L) /S*CMPLX ( COS ( ARG ) ,SIN(ARG) ) 
UB=XT ( L) *UA 
DO  25  M=M3 , M4 
FA (M) =BJ (M) *UA+FA (M) 

FB(M)=BJ(M) *UB+FB (M) 

25  CONTINUE 
15  CONTINUE 

IF(Ml.NE.O)  GO  TO  26 
FA ( 1 ) =FA ( 3 ) 

FB ( 1 ) =FB ( 3 ) 

26  UA=U1 

DO  27  M=M5 , M6 

M7=M-1 

M8=M+1 

F2A=UA* ( FA (M8 ) +FA (M7 ) ) 

F2B=UA* (FB(M8) +FB(M7) ) 

UB=U*UA 

F1A=UB* (FA (M8) -FA (M7) ) 

F1B=UB* (FB(M8) -FB(M7) ) 

U4=D4  *UA 

U2=D3*F1A+U4*FA(M) 

U3=D3*F1B+U4*FB(M) 

U4=DR*F2 A 

U5=DR*F2B 

K1=K2-1 

K4=K1+N 

K5=K2+N 

R (K2+MT) =-D5* (F2A+D2*F2B) 

R (K5+MT) =D1* (F1A+D2*F1B) 

IF(IP.EQ.l)  GO  TO  21 
R(K1) =R(K1) +U2-U3 
R (K4 ) =R (K4 ) +U4-U5 
IF(IP.EQ.MP)  GO  TO  22 
21  R(K2)=U2+U3 
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5 

6 

7 

8 
9 

10 

11 

12 

13 

14 

15 

16 

17 

18 

19 

20 
21 
22 

23 

24 

25 

26 

27 

28 

29 

30 

31 

32 

33 

34 

35 

36 

37 

38 

39 

40 

41 

42 

43 

44 

45 

46 

47 

48 

49 

50 

51 


R(K5)=U4+U5 
22  K2=K2+N2 
UA=UB 

27  CONTINUE 

12  CONTINUE 
R(NPl-l) =(0. ,  0.  ) 

R (NP1 ) = ( 0 . ,0. ) 

R (NP1+MT)  =  ( 0 .  ,  0.  ) 

RETURN 

END 

SUBROUTINE  PLAN (Ml , M2 , NP , NT , RH , ZH , XT , AT , THR , R) 

C 

C  PLANE  WAVE  RECEIVE  VECTOR  ELEMENTS.  CAN  BE  USED  TO  COMPUTE 
C  THE  PLANE  WAVE  EXCITATION  ELEMENTS  IF  PLANE  WAVE  SCATTERING 
C  IS  BEING  CALCULATED.  (THE  RELATIONSHIP  BETWEEN  THE  RECEIVE 
C  AND  EXCITAION  ELEMENTS  AS  DISCUSSED  BY  MAUTZ  AND 
C  HARRINGTON.) 

C 

COMPLEX  R ( 1800 ) , U , U1 , UA, UB , FA ( 5 ) ,FB(5) , F2A, F2B, F1A, FIB 
COMPLEX  U2 ,  U3  , U4 ,  U5 , CMPLX 

DIMENSION  RH(1000) , ZH(1000) , XT ( 10) , AT ( 10) , R2 ( 10) 
DIMENSION  Z2 (10) , BJ (5000) 

MP=NP-1 
MT=MP-1 
N=MT+MP 
N2  =  2  *N 
CC=COS (THR) 

SS=SIN (THR) 

U=(0. , 1 . ) 

Ul=3. 141593*U**M1 

M3=M1+1 

M4=M2+3 

IF(Ml.EQ.O)  M3  =  2 

M5=Ml+2 

M6=M2+2 

DO  12  IP=1 , MP 

K2  =  IP 

I=IP+1 

DR= . 5  * ( RH ( I ) -RH (IP) ) 

DZ=. 5* ( ZH ( I ) -ZH(IP) ) 

D1=SQRT ( DR*DR+DZ*DZ ) 

Rl= . 2  5* (RH ( I ) +RH (IP) ) 

IF (R1 . EQ. 0 . )  Rl=l . 

Z1=.5*(ZH(I)+ZH(IP) ) 

DR= . 5* DR 

D2=DR/R1 

DO  13  L= 1 , NT 

R2 ( L) =R1+DR*XT ( L) 

Z2 (L)=Z1+DZ*XT(L) 

13  CONTINUE 
D3=DR*CC 
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1 

D4=-DZ*SS 

2 

D5=D1*CC 

3 

DO  23  M=M3,M4 

4 

FA (M) =0 . 

5 

FB (M) =0 . 

6 

23 

CONTINUE 

7 

DO  15  L=1 , NT 

8 

X=SS*R2 (L) 

9 

IF (X . GT . . 5E-7 )  GO  TO  19 

10 

DO  20  M=M3,M4 

11 

BJ (M) =0 . 

12 

20 

CONTINUE 

13 

BJ  ( 2 )  =1 . 

14 

S=l. 

15 

GO  TO  18 

16 

19 

M=2 . 8*X+14 . -2 ./X 

17 

I F ( X . LT . . 5 )  M=ll . 8+ALOG10 (X) 

18 

IF(M.LT.M4)  M=M4 

19 

BJ (M) =0 . 

20 

JM=M- 1 

21 

BJ ( JM) =1 . 

22 

DO  16  J=4 , M 

23 

J2=JM 

24 

JM=JM- 1 

25 

J1=JM-1 

26 

BJ ( JM) =J1/X*BJ ( J  2 ) -BJ ( JM+2 ) 

27 

16 

CONTINUE 

28 

S  =  0. 

29 

IF (M . LE . 4 )  GO  TO  24 

30 

DO  17  J=4,M,2 

31 

S=S+BJ ( J ) 

32 

17 

CONTINUE 

33 

24 

S=BJ ( 2 ) +2 . *S 

34 

18 

ARG=Z2 (L) *CC 

35 

UA=AT (L) /S*CMPLX ( COS ( ARG) , SIN(ARG) ) 

36 

UB=XT ( L) *UA 

37 

DO  25  M=M3 , M4 

38 

FA (M) =BJ (M) *UA+FA (M) 

39 

FB (M) =BJ (M) *UB+FB(M) 

40 

25 

CONTINUE 

41 

15 

CONTINUE 

42 

IF(Ml.NE.O)  GO  TO  26 

43 

FA ( 1 ) =FA ( 3 ) 

44 

FB ( 1 ) =FB ( 3 ) 

45 

26 

UA=U1 

46 

DO  27  M=M5,M6 

47 

M7=M-1 

48 

M8=M+1 

49 

F2A=UA* (FA (M8 ) +FA (M7 ) ) 

50 

F2  B=UA* ( FB (M8 ) +FB (M7 ) ) 

51 

UB=U*UA 
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1 

F1A=UB* ( FA (M8 ) -FA (M7 ) ) 

2 

F1B=UB* (FB(M8) -FB(M7) ) 

3 

U4  =  D4  *UA 

4 

U2=D3*F1A+U4*FA(M) 

5 

U3=D3*F1B+U4*FB(M) 

6 

U4=DR*F2A 

7 

U5=DR*F2  B 

8 

K1=K2-1 

9 

K4=K1+N 

10 

K5=K2+N 

11 

R (K2+MT) =-D5* (F2A+D2*F2B) 

12 

R (K5+MT) =D1* (F1A+D2*F1B) 

13 

IF ( IP . EQ . 1 )  GO  TO  21 

14 

R ( K1 ) =R ( K1 ) +U2-U3 

15 

R (K4 ) =R ( K4 ) +U4 -U5 

16 

IF ( IP. EQ . MP)  GO  TO  22 

17 

21 

R(K2) =U2+U3 

18 

R (K5) =U4+U5 

19 

22 

K2=K2+N2 

20 

UA=UB 

21 

27 

CONTINUE 

22 

12 

CONTINUE 

23 

RETURN 

24 

25 

END 
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